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Abstract 


Three mechanisms which can induce thermo-hydraulic oscillations at 
near-critical and at super-critical pressures are distinguished and dis- 
cussed. 

Experiments show that low frequency flow oscillations are most pre- 
valent in systems of practical interest. A quantitative formulation and 
analysis is therefore presented concerned with predicting the onset of 
these "chugging" oscillations as function of fluid properties , system 
geometry and operating conditions. 

The problem is analyzed by perturbing the inlet flow, linearizing 
the set of governing equations and integrating them along the heated duct 
to obtain the characteristic equation. The latter is given by a third 
order exponential polynomial with two time delays. 

Conditions leading to aperiodic as well as to periodic flow 
phenonema are investigated. The first pertains to the possibility of 
flow excursion the latter to the onset of flow oscillation . 

Stability maps and stability criteria are presented which, previously, 
were not available in the literature. They can be used to determine: 

a) The region of stable and unstable operation and 

b) The effects which various parameters have on either promoting 
or preventing the appearance of flow oscillations. 

The effects of various parameters are analyzed and improvements are 
suggested whereby the onset of flow oscillation can be eliminated. 

The similarity between "chugging" combustion instabilities and 
thermally induced flow oscillations at near and super-critical pressures 
is pointed out. 

A review of the present understanding of the near-critical thermo- 
dynamic region is also presented. 
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1 . Research Objectives 


This research was conducted to determine the fundamental nature of 
oscillation, and of instabilities in the flow of cryogenic fluids with heat 
addition. 

The investigation was motivated by the fact that severe oscillations 
have been experienced in rocket engines heat exchangers utilizing oxygen 
and hydrogen at both subcritical and supercritical pressures. 

The particular objectives of this investigation weres 

1. To distinguish a number of mechanisms which may be respon- 
sible for thermally induced flow oscillations at near cri- 
tical and at supercritical pressures. 

2.. To present a quantitative formulation of the mechanisms 

which appear, to be most significant from the point of sys- 
tem design and operation. 

3. To predict the onset of these oscillations in terms of the 
geometry and of the operating condition of the system. 

4. To analyze the consequences of the theoretical predictions 
and to suggest improvements whereby the onset of these 
oscillations can be eliminated. 


2. Summary and Conclusions 

1. Mechanisms Leading to Unstable Operation 

Three mechanisms which can induce thermo -hydraulic oscillations at near 
critical and at supercritical pressures have been distinguished. 

One is caused by the variation of the heat transfer coefficient at the 
transposed i.e., at the pseudo-critical point. 

The second is caused by the effects of large compressibility in the 
critical thermodynamics region. 

Finally, the third mechanism is caused by variations of flow character- 
istics brought about by variations of fluid density during the heating pro- 
cess. The propagations of these variations through the system introduce 
various time delays which, under certain conditions, can cause unstable flow 

This last mechanism, which induces low frequency oscillations, was 
investigated in detail because available experimental data show that this 
type of flow oscillations is most prevalent in systems of practical interest 

2 o Formulation of the Problem 

The problem was formulated in terms of an equation of state and of 
three field equations describing the conservation of mass, energy and mo- 
mentum. 

The subcritical pressure range of operation was differentiated from the 
supercritical one by using the appropriate equation of state. 

The problem was analyzed by perturbing the inlet flow, linearizing the 
set of governing equations and integrating them along the heated duct to 
obtai n the characteristic equation. 



• The Characteristic Equation 



The characteristic equation is given by a third order exponential- poly- 
nomial with two time constants, (see Eq. V-15). It is expressed in terms of 
fluid properties, of system geometry and of operating conditions by means of 
influence coefficients (see»Eq. V-16 through Eq. V-22). 

The influence coefficients express the effects of the inlet flow 
perturbation and of the space lag perturbation on the various pressure drops 
of the system. By introducing various definitions for the average, for the 
log^mean and for the mean densities and velocities, it is shown that each 
pressure drop is weighed with respect to a different velocity . This 
result, which follows, from the integration of the governing set of 
equations, i*e., from the distributed parameter analysis, could not have 
been obtained from an analysis, based on "lumped" parameters. Consequently 
the accuracy of an analysis based on this latter approach can be estimated 
by means of the results obtained in this investigation* 

The characteristic equation was used to obtain stability maps and 
stability criteria which, previously, were not available in the literature. 
The stability maps and criteria can be used to determine 

a. The region of stable and of unstable operation and 

b. The effects which various parameters may have on either promoting 
or on preventing the appearance of flow oscillations. 

Conditions leading to aperiodic as well as to periodic flow phenomena 
were investigated. The first pertain to the possibility of flow excursion 
whereas the second pertain to the onset of flow oscillation . For this latter 
case the flow stability . in systems with low inlet subcodlihg was considered 
separately from that corresponding to systems with high iiilet subcooling- 

III 


The stability problem at intermediate subcooling will be considered in a 
future report- 

4» Excursive Flow Instability 

It was shown that, at supercritical pressures, a flow system with heat 
addition can undergo flow excursions because the hydraulic characteristics 
of the system are given by a cubic relation between the pressure and the 
mass flow rate (see Eq° VI-20)- The latter is a consequence of density 
variations in the system. 

This excursive flow instability, at supercritical pressures, is the 
equivalent of the "Ledinegg" excursive instability in boiling systems at 
subcritical pressures- This equivalence is supported by experimental data 
(see Figure VI-1) which show : ! that in both pressure regions, the flow system 
has similar hydraulic characteristics. 

A stability criterion which predicts the onset of the excursive in- 
stability was derived in terms of system geometry, of fluid properties and 
of operating conditions, i-e-, of system pressure, flow rate, inlet temp- 
erature and power input (see Eq- VI-13) ° Various aspects of this type of 
instability are discussed together with provisions required to prevent its 
appearance (see Section VI-2). 

5- Flow Oscillations at Low Inlet Subcooling 

It was shown that for a system with low inlet: subcooling the character- 
istic equation can be reduced to a second order polynomial with one time 
delay (see Eq° VII-7). For such a system the propensity to flow oscillation 
can be analyzed by means of the stability maps recently presented by Bhatt 
and H su (see Figure VII- 1.) » ' 


It was shown further that, when the inertia can be neglected in a system 
with low inlet subcooling then the characteristic equation reduces to a first 
order exponential polynomial with one time delay (see Eq. VII-16). 

For such a system the flow will be unconditionally stable if the 
stability number N s (defined by Eq. VII-39) is larger than unity. If 
N s is smaller than unity, stable operation is still possible ij: the angular 
frequency of the inlet perturbation is larger than the critical one (given 
in Eq. VII-40) or if the transit time is shorter (than the critical one 
(given by Eq. VII-41) . 

The region of stable and of unstable operation are shown in a stability 
map (see figure VII-2) which can be used to analyze the effects that various 
parameters have on the propensity to induce or to prevent flow oscillations 
(see Section VII. 3). 

Although the analytical predictions have not yet been tested quantita- 
tively, the trends predicted by this map and by the stability criterion 
(see Eq. VII-22 or Eq. VII-29) are in qualitative agreement with experimental 
observations (see Section VII. 3). 

6. Flow Oscillation at High Inlet Subcooling 

It was shown that when the effects of the two time delays can be 
neglected then the characteristic equation reduced to a * thir'O ..ordex 1 polynomial 
(see Eq. VIII-2). A stability criterion was also derived (see Eq. VIII-15) 

* r 

which specifies the conditions for stable operation. 

Various aspects of this type of oscillations were discussed together 
with provisions required to prevent their appearance (see Section VIII-2). 

It is shown that the flow is more stable at high subcoolings than at low. 
Furthermore, it is concluded that the destabilizing effect of subcooling 


must go through a maximum at intermediate range, (see Section VIII-2). 

7 » Significance of the Results 

The results of this analysis indicate several improvements in the design 
and/or in the operating conditions which can be made to prevent the onset of 
flow excursions or of flow oscillations. These are discussed in more detail 
in relation to each type of instability (see Sections VI-2, VII-3, and 
VIII-2). 

It was shown that the predominance of a particular parameter results 
in a particular wave form and in particular frequency (see Eq. VXI-40 and 
Eq. VIII-17). This result indicates that the primary cause of the instability 
can be determined from the trace of flow oscillations. 

Perhaps the result of greatest significance revealed in the present 
investigation is the similarity between the characteristic equation which 
predicts ''chugging" combustion instabilities and the characteristic equation 
which predicts the thermally induced flow oscillations for fluids in the 
near critical and in the supercritical thermodynamic region. Since it is 
well known that "chugging" combustion instabilities can be stabilized by an 
appropriate servo-control mechanism, the results of this investigation 
indicate that low frequency flow oscillations, at near critical and at 
supercritical pressures may be also stabilized. 

The preceding conclusions have not yet been tested against experimental 

* 

data. If confirmed, then the results of this study will provide a method 
whereby stable operation can be insured in an in trinsicalLy unstable region. 


3. Recommendations 


The recommendations listed in the four tasks below^ define the effort 
needed to complete and to verify the results obtained in this investigation. 

1. Verify the stability criteria based on the second and first order 
exponential polynomials which have been derived in the course of 
these investigations. For this purpose use available experimental 
data for various fluids at subcritical and at supercritical pres- 
sures. 

2c From the characteristic equation given by the third order exponential, 
polynomial with two time delays (Eq. V-15) derive stability maps and 
stability criteria applicable to the entire range of subcoolings. 

Test these results against available experimental data. 

3. Modify the characteristic equation to take into account fhe effects 
of the entire flow system i.e. ; of the flow loop. In particular in- 
clude the effects of the inertia of the liquid in the storage tank and 
in the supply lines together with the flow and elastic characteristics 
of these lines. 

4. B.ased on the results obtained from the preceding three tasks specify a 
servo-control mechanism which could be used to stabilize the flow for 
a system of practical interest and verify the results by experiments* 
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I. Introduction 


I . 1 The Problem and Its Significance 

A fluid in the vicinity of the critical point is an efficient heat 
transfer medium because of the large specific heat and of the large co- 
efficient of thermal expansion,, Consequently;, the demand for increased 
efficiency of several advanced systems generated an interest in employing 
fluids at critical and supercritical pressures either as cooling or working 
media . For example , nuclear rockets, power reactors, high pressure once- 
through boilers, regenerative heat exchangers for rocket engines and a new 
sea water desalinization process are designed to operate in the critical 
and the supercritical thermodynamic region. These developments made it 
necessary to obtain data on and to improve the understanding of the thermal 
and the flow behavior over a broad range of fluid states, 

A great number of investigations conducted for such a purpose have 
revealed that, in the critical as well as in the supercritical thermo- 
dynamic region, flow and pressure oscillations may occur when certain 
operating conditions are reached. These oscillations were observed in 
systems with forced flow as well as with natural circulation. 

The occurrence of sustained pressure and flow oscillations and the 
attendant temperature oscillations are very undesirable and detrimental to 
reliable operation of a system. Mechanical vibrations and thermal fatigue 
induced by these oscillations very often result in a rupture of the duct. 

In liquid propellant rocket engines flow and pressure oscillations can also 

induce combustion instabilities resulting in a breakdown of the system. 
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Furthermore, in nuclear reactor systems flow and pressure oscillations may 
induce divergent power oscillations leading to the destruction of the 
entire system. Consequently, there is considerable practical interest and 
incentive to investigate, quantitatively, the conditions leading to the 
inception of these oscillations, 


I. 2 Previous Work 

Severe pressure and flow oscillations were observed in experiments 
performed with various fluid in the supercritical thermodynamic region. 

Such oscillations were reported by Schmidt, Eckert andGrigull til. (ammonia) 
Goldman j^2, 3^ , (water); Firstenberg ^4^ , (water); Harden j~5~| , (Freon-114); 
Harden and Boggs , (Freon-114); Walker and Harden £ 7~] , (water, Freon-114, 
Freon-12, carbon dioxide); Holman and Boggs [_8~j, (Freon-12); Hines and 
Wolf |^9 1 (RP-1 and diethycyclohexane) ; Platt and Wood J"lO~j (dxygen); 
Ellenbrook, Livingood and Straight jjLl^, (hydrogen); Thurston £l2~j, (hydrogen, 
nitrogen); Shitzman jjL3, 14*^ (water); Semenkover jjL5^ (water) ; Cornelius and 
Parker (Freon-114); Cornelius ^17~| (Freon-114); and Krasiakova and 

Glusker j^l8"j (water). 

For a given fluid the characteristics (frequency and amplitude) of 
these oscillations varied with operating conditions. In general, two types 
of oscillations were observed: acoustical and chugging oscillations. For 

example, Shitzman £l3~j reports that, for water at 250 atm, the pressure and 
temperature oscillations had a period of 80 sec. and an amplitude of 25 atm. , 
and of 100°C respectively. Decreasing the flow rate and the power density 
resulted in decreasing the period to 15 sec. However, at a pressure of 
5000 psi, Goldman J^2, 3 reports pressure oscillations with frequencies 
from 1400 to 2200 cps. Similar high frequency oscillations (1000-10,000 cps 


and 380 psi peak to peak) were reported by Hines and Wolf £ 9*1 for RP-1, 

Three classes of pressure oscillations in the supercritical region 
were observed in the experiments of Thurston jj^j; these were described ass 

1) Open-open pipe resonance observed at medium and high flow rate,. 

This mode is associated with the fundamental wavelength of an 
open-open pipe„ 

2) Helmholtz resonance g associated with a resonator composed Of a 
cavity connected to an external atmosphere via an orifice or neck, 

3) Supercritical oscillations appearing usually at low flow rates, 

Hines and Wolf £ 9*"J, however , report only two general types of oscilla- 
tions o a high frequency (3000-75000 cps) oscillations audible as a clear 
and steady scream and an oscillation with a lower frequency (600-2400 cps) 
which was audible as a chugging or pulsating noise. The dominant frequencies 
of these oscillations did not correspond to simple acoustic resonant fre- 
quencies for the tubes. 

Cornelius and Parker^l6, 17~j describe in detail the two types of 
oscillations and note that the frequency of the acoustical oscillations 
decreases with temperature whereas the frequency of the chdgging oscillations 
increases with temperature. Occasionally, both types of oscillations occured 
simultaneously, 

A quantitative formulation and explanation of the conditions leading 
to the. appearance of the pressure and flow oscillations has not been reported 
yet, although several qualitative explanations have been advanced. It is 
generally agreed that the oscillations are caused by the large variations 
of the thermodynamic and the transport properties of the fluid as it passes 
through the critical thermodynamic region. 


Several investigators (12, 13, 14) note that the appearance of oscil- 


lations occurs when the temperature of the heating surface exceeds the 
" pseudo critical " or the " transposed " critical temperatures, i.e., the tem- 
perature where the specific heat reaches its maximum value (see Figure A1 in 
Appendix A). Oscillations were not observed if the inlet temperature was 
above this temperature. From this it was concluded that the mechanism for 
driving the oscillation occurred only when a "pseudo liquid" state was present, 
in some parts of the heated duct. 

Firstenberg (4) attributes the oscillations to the variations of the heat 
transfer rates to the fluid, whereas Goldman (2, 3) explains the oscillations 
as well as the steady state heat transfer mechanism in the critical and super- 
critical thermodynamic region as "boiling like" phenomena associated with non- 
equilibrium conditions. According to Goldman, below the pseudo-critical temper- 
ature the fluid is essentially a liquid, above this temperature it behaves as a 
gas. At the pseudocritical temperature, the density gradient and the specific 
heat reach maximum values giving an indication of the energy required to over- 
come the mutual attraction between the molecules. The fluid in the immediate 
vicinity of the heated wall is in a gas -like state; whereas the bulk fluid may 
still be in the liquid-like state. If by means of turbulent fluctuations a 
liquid-like cluster is brought into contact with the heating surface a large a- 
mount of energy will flow from the surface to the cluster because of the -large 
temperature difference and because of the high conductivity of the liquid-like 
fluid. This energy may be large enough to "explode" clusters of molecules from 
the liquid-like state to the gas -like state. Thus, according to Goldman (2, 3), 
one may visualize the supercritical region as a region where explosions of liquid- 
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like clusters into gas-like aggregates take place , Goldman considers this 
process to be similar to the formation of bubbles in liquids during boiling 
at suberifcical pressures , 


The conditions under which oscillations occur were summarized by- 
Goldman as follows; 

1) Heat transfer with "'whistle" (ioe. s with oscillations) occurs 
only at high heat flux densities and with bulk temperatures lower 
than the pseudocrttical temperature, 

2) At a given flow rate and inlet temperatures whistles occur at 
higher flux densities for higher pressure levels, 

3) At given flow rate and pressure^ whistles occur at lower heat 
flux densities for higher inlet temperatures, 

4) At a given pressure and inlet temperature 9 whistles occur at 
higher heat flux densities for higher flow rates, 

5) Whistles can be produced with various lengths of the test sections 
but the heat flux or inlet temperature must be increased to bring 
it about if the tube is shortened, 

Visual observation that boiling-like phenomena can exist at supercritical 
pressures was reported by Griffith and Saberski [l9] in e&£>eriments conducted with 
It- 114 , The photographs of the process revealed a behavior similar to 
that observed in pool boiling at subcritical pressures. 

Similarly 5 high speed movies of hydrogen at supercritical pressures 
taken by Graham s et al £ 20 *"| revealed a phenomenon resembling boiling. Clusters 
of low density fluid were observed rising through a denser fluidrgiving 
boiling- like -appearance^- 



Hines and Wolf ^9*( attribute the appearance of the flow oscillations 

4 

at supercritical pressures to the variations of liquid viscosity,, They 
note that a small change of temperature near the critical point results 
in a large change -of viscosity. Consequently, a sudden increase in wall 
temperature could cause a thinning of the laminar boundary layer due to 

variation of the viscosity. Thinning of the boundary layer would result 

• * *•*' .* 

in a drop of the wall temperature and a corresponding increase of viscosity. 
This would cause a thicker boundary layer and produce another rise of tem- 
perature, thus repeating the cycle. It was shown by Bussard and DeLauer [j2lj 
and by Harry £ 22^ that a viscosity-dependent mechanism can induce an unstable 
flow in single phase flow systems when the absolute gas temperature .is in- 
creased by a factor of 3.6 or more. ’ Such flow oscillations were observed 
by Guevara et al^23~|with helium flowing through a uniformly heated 
channel. 

Harden and co-workers ^5 , 6, 7*] concluded from their experiments that 
sustained pressure and flow oscillations appeared when the bulk fluid reached 
a temperature at which the product of the density and enthalpy has its maximum 
value. This explanation was, however, criticized by Cornelius ^17.]. 

Cornelius and Parker |^16, 17~j postulate that both acoustical and the 
chugging oscillations originate in the heated boundary layer. When the 
fluid in the heated boundary layer is in a "pseudo vapor" state, a pressure 
wave passing the heated surface would tend to compress the boundary layer, 
improve the thermal conductivity and cause an increase of the heat transfer 
coefficient. A rarefraction wave passing over the heated wall would have 
just the opposite effect. Thus, this pressure--dependent — heat-transfer 
rate could induce and maintain an acoustic oscillation. Cornelius and Parker 
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attribute the appearance of chugging oscillations to "boiling -like" phenomena 
and a sudden improvement of the heat transfer coefficient. An approximate 
numerical solution verified the importance of the heat transfer improvement 
in triggering and maintaining oscillations. 

Of particular interest to the analysis presented in this paper are the 
experimental results of Semenkover, (JL53 and of Krasiakova and Glusker £ 18^ « 
for water at 250 atm. For a constant power input Q to the system their data 
show a pressure versus mass flow relation that is illustrated in Figure 1-1. 

It can be seen that for large values of inlet enthalpy i, there is a monotonic 
increase of pressure drop with flow rate. At a certain lower value of i, 
the curve shows an inflection point. For still lower values of inlet enthalpy, 
there is a region where the pressure drop decreases with increasing flow rate. 

Such a pressure drop=flow rate relation occurs in boiling systems and gives 
rise to an excursive type of instability which was analyzed first by Ledinegg 
L 24 J and by numerous investigators since £25 - 4?J 0 Consequently, the data 
of {l5, 18 J tend to confirm the similarity between instabilities. observed 
during subcritical boiling and those observed at supercritical pressure 
suggesting therefore a common mechanism. 

1.3 Purpose and Outline of the Analysis 

From the preceeding brief review of the present understanding of flow 
oscillations at supercritical pressures, it can be concluded that several modes of 
oscillation exist. If can be expected, therefore, that several mechanisms can 
be effective in inducing unstable flow. Indeed, as discussed in the preceeding 
section, several qualitative explanations of the phenomenon have been already 
advanced. However, a quantitative formulation of the problem is still lacking. 



I 




YDRAULIC CHARACTERISTICS FOR WATER AT SUPERCRITICAL 
RESSURE ( P = 250 atm) FLOWING THROUGH A HEATED DUCT. 
ATA OF SEMENKOVER FOR VARIOUS INLET ENTHALPIES: 

- 400, 2- 350,3-300,4- 200,5-100 IN kcal/kg. 





The analysis presented in this paper has four objectives; 

1) To distinguish a number of mechanisms which may be responsible for 
thermally induced flow oscillations at nearcritical and at supercritical 
pressures 

2) To present a quantitative formulation of the mechanism which appears to 

be most significant from the point of view of system design and operation. 

3) To predict the onset of these oscillations in terms of the geometry and 
of the operating conditions of the system. 

4) To analyze the consequences of the theoretical predictions and to suggest 
improvements whereby the onset of these oscillations can be eliminated. 

The particular mechanism which is formulated and analyzed in this paper 
is based on the effects of time lag and of density variations. It is well 
known that these effects can induce combustion instabilities in liquid pro- 
pellant rocket motors as discussed by Crocco and Cheng 01 - It was shown 
by Profos^49^ , Wallis and Heasley ^ 50)"”^ and by Bourd £si"^ that the effects 
of time lag and of density variations can also induce unstable flow in 
boiling mixtures at subcritical pressures. The suggested similarity of 
flow oscillations observed at supercritical pressures with those observed in 
two phase mixtures at subcritical pressures prompts us to formulate and 
to analyze the problem in terms of this mechanism. In particular, the ex- 
perimental results of Semenkover ^15*) and of Krasiakove and Glusker 18^ 
discussed in the preceding section., together with the chugging oscillations 
described by several authors provide enough evidence to warrant a more de- 
tailed analysis of flow oscillations at supercritical pressures in terms 
of the time lag effect. 

The present analysis is similar to those reported by Wallis and Heasley 


rsoj and by Boure [5lJ in two respects : the formulation and the assumptions 

are the same. In particular s it is assumed that the density of the medium 
is a function of enthalpy only. The effects of pressure variations are s 
therefore neglected.* As noted by Wallis and Heasley D<0 this assumption 
results in the decoupling of the momentum equation from the energy and the 
continuity equations. The momentum equation can be integrated then separately 
after the velocity and density variations are obtained from the continuity and 
the energy equations. Following Bour4 £ 51 J the problem is formulated in terms 
of an equation of state and of three field equations describing the conser^ 
vat ion of mass s energy and mementum. 

Apart from the fact that the analyses of Wallis and Heasley CsoJ and 
of Bounl [5lJ were derived to predict unstable flow in boiling two phase 
mixtures the present analyses (concerned with flow oscillations at near= 
critical and at supercritical pressures) differs from t Heirs in two respects? 

1) the form of the constitutive equation of state is different „ 2) the 
characteristic equation describing the onset of oscillations is different. 

From this characteristic equations we shall derive stability maps and 
stability inertia which s previously, were not available in the literature. 

The outline of the paper is as follows. In Chapter II some general 
comments are made regarding 1) the nature of the thermally induced flow 
oscillations at nearcritical and at supercritical pressures^ 2) tne effect 
of tne time lag 8 3) the implication and limitations of tne assumptions and 
4) the general method of solution. In Chapters III and IV the problem is 
formulated and the set of governing equation is solved. The characteristic 
equation which predicts tne onset of flow oscillation is derived in Chapter Vj 
it is of tne form of a third order exponential polynomial with two time 
delays. From the characteristic equation a criterion is derived in 

*The limitations and implications of this assumption are discussed in Chapter II 
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Chapter VI which predicts the onset of an excursive type of instability at 
supercritical pressures.* This excursive instability at supercritical 
pressure is the equivalent of the so called Ledinegg excursive instability 
for boiling at subcritical pressures. The effect of time lags in inducing 
flow oscillations is analyzed in Chapters VII and VIII which consider first 
and second order expotential polynomials. Stability diagrams which predict 
the regions of stable and unstable flow in terms of the operating parameters 
are given in these two chapters together with suggested improvements whereby 
the onset of ..oscillations can be eliminated. The recommendations for 
future work and the conclusions are given in Chapters IX . 

The status of the present understanding of thermodynamic phenomena that take 
place in the critical thermodynamic region is discussed in Appendix A. 

1.4 The Significance of the Results 

Three mechanisms which can induce thermo -hydraulic oscillations at 
supercritical pressures have been distinguished in this paper. One is 
caused by the variation of the heat transfer coefficient at the transposed s 
i.e.j, pseudo critical point. The second is caused by the effects of large 
compressibility and the resultant low velocity of sound in the critical 
region. Finally, the third mechanism is caused by the large variation of flow 
characteristics brought about by density variations of the fluid during the 
heating process. The propagations of these variations in particular of the 
enthalpy and of the density, through the system introduce delays which. 


*This criterion was first derived by the writer in the Second Quarterly 
Progress Report, "Investigation of the Nature of Cryogenic Fluid Flow 
Instabilities in Heat Exchangers," Contract NAS8-11422, 1 February 1965. 
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under certain conditions, can cause unstable flow* This last mechanism, that 
induces low frequency oscillations, is investigated in detail because ex- 
perimental data show that this type of oscillation is most prevalent . 

It is shown that at supercritical pressure unsteady flow qonditions 
both excursive and oscillatory can occur. A characteristic equation is 
derived that predicts the onset of flow instabilities caused by density 
variations in the critical and supercritical thermodynamic region. The 
same characteristic equatipn can be used to predict the onset of flow 
instabilities in boiling at subcritical pressure, jlf the effect of the 
relative velocity between the two phases can be neglected. Experimental 

evidence shows that this effect becomes negligible at reduced pressures 

/ 

above say 0.85. Consequently, at ne ar critical and supercritical pressures, 
the characteristic equation, which is expressed in terms of system geometry 

f. 

and operative conditions, can be used to determines 

a) The region of stable and unstable behavior. 

b) The effect which various parameters may have on either promoting 
or on preventing the appearance of flow oscillations. 

From this characteristic equation simple "rule of thumbs" criteria are 
also derived based on the assumption that one or the other of the various 
parameters is dominant. It is shown that the dominance of a particular 
parameter results in a particular frequency and wave form. This results 
permits a diagnosis of the primary cause of the instability from the trace 
of flow oscillations. 

It is of particular interest to note that- the characteristic equation 
derived in this paper for predicting flow oscillations at supercritical 
pressure is of the same form as a characteristic equation derived by Crocco 



and ChengC48]to predict combustion instabilities of liquid propellant 
rocket motors,,* It is well established in the combustion literature that 
a servo-control mechanism can be used to stabilize the low frequency com- 
bustion instability. The similarity of the characteristic equations is 8 
therefore 9 significant because it indicates that stable operation could 
be insured also in the nearcritical and in the supercritical region by 
using an appropriately designed servo-control mechanism. 


*This similarity between combustion and two phase flow instabilities should 
not come as a surprise if one recalls that the processes of combustion and 
of boiling are both chemical processes involving large enthalpy and density 
changes. 
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II o General Considerations 


II o 1 The System and the Thermodynamic Process 

In order to understand the mechanisms of the thermally induced flow 
oscillations at supercritical pressures s it is necessary to examine 
briefly the system and the thermodynamic process . 

The system of interest is shown in Figure 11=1. It consists, of a 
fluid flowing through a heated duct of length 1, Without loss of 
generality it will be assumed that the duct is uniformly heated at a 
rate of Q. Two flow restrictions are located one at the entrances the 
other at the exit of the duct. 

The thermodynamic process starts with the fluid at £ supercritical 
pressure P 3 entering the heated duct with velocity U t , The temperature 
Tp of the fluid at the inlqt is well below the critical temperature of the 
fluid under considerations. As the energy is being transferred from the 
heated duct to the fluid its temperature T 3 specific volume v 3 and enthalpy 
i, will increase. Thus., the temperature T^, at the exit may be considerably 
above the critical temperature. In a number of systems of practical interest 
it can be assumed that this process takes place at an approximately constant 
pressure. 

In order to formulate the problem it is necessary to specify the 
constitutive equation of state which describes the relation between say 
the. specific volume, the pressure and the enthalpy for the particular 
fluid. This requires data oh the thermodynamic properties of the fluid in 
the region of interest. The region of interest to this study are the 
nearcrftical and the supercritical regions. 




The present understanding of the thermodynamic properties and phenomena 

at nearcritical and supercritical pressures is. reviewed in Appendix A. 

It is shown there that at these pressures the fluid has the characteristics 

of a liquid when the temperature is sufficiently below the critical one. 

However, if the temperature is increased sufficiently above the critical 

temperatures the fluid will have the characteristics of a gas. This is 

illustrated in Figure II-2 which is a plot of the specific volume and of 

the temperature versus the enthalpy for. oxygen at a reduced pressure of 

P -1.1. 

r 

It can be seen from this figure that at low enthalpies the specific 
volume is essentially constants this is a characteristic of liquids. As 
the enthalpy increases the specific volume increases approaching values 
predicted by the perfect gas law. It can be seen also that this change 
from a liquid- like state (region © - ©) to a gas-like state (region 
0 - ©) occurs over a transition region denoted by © - © 
on Figure II-2. 

It appears, therefore 9 that at supercritical pressures the relation 
between the specific volume and the enthalpy can be approximated by con- 
sidering three regions; a liquid-like , a transition and a gas like region. 
For oxygen Figure II-2 indicates also that., .as a first approximation, the 
transition region can be reduced to a transition point, reducing, therefore, 
the problem to a "two-region" approximation.* Since oxygen is the fluid 
of primary interest to this analysis, we shall use the "two-region" 

*The "three region" approximation was first introduced by the writer in 
analyzing the excursive instability at supercritical pressures (see foot- 
note on page 10). Following this work Dr. R. Fleming, from the Research 
and Development Center of the 6E Co. , introduced the "two region" approxi- 
mation for oxygen. These two approximations are discussed further in 
Appendix B. 
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approximation for describing the relation between the specific volume 
and the enthalpy, It is assumed , therefore s that the l? heavy' ? fluid (of 
constant density) persists until the transition point is reached; above 
this point the fluid will have the properties of the “light 8 phase. It 
remains now to define this transition point. 

In boiling at subcritical pressure the transition from the heavy to 
the light phase corresponds to the onset of boiling. Consequently 9 it 
will be assumed that in the nearcritical region the transition point 
corresponds to the enthalpy at saturation temperature. 

At supercritical pressures it will be assumed that the transition point 
corresponds to the transposed, critical points i,e.» to the pseudo critical 
point which is defined as the point where C^ reaches its maximum value. 

It is discussed in Appendix A that the locus of pseudo critical points 
can be regarded as the extension of the saturation line in the super- 
critical region. 

IT. . 2 Time Lag and Space Lag 

It is of interest to consider now the timewise and spacewise des- 
cription of the process.* 

If we follow a particle from the time it enters the heated section 
until it leaves it, we shall observe that its properties change from v^ 

— . - | — - — A 

*We .follow here Crocco and Cheng who gave an equivalent description 

for combustion instabilities. The same comment applies to the three 
sections that follow. Indeed, this reference proved to be most stimu- 
lating and useful in the course of this investigation. 
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and i^ at the inlet to v^ and i 3 at the exit (see Figure 11=3). In 
view of the "two=region" approximation we would note that the transition 
from the "heavy'* to the "light” fluid occurred when the properties (specifically 
the enthalpy) reached values that correspond to the transition point. The 
time elapsed between the injection of the heavy particle in the heated duct 
and its transformation to the "light" fluid will be denoted as the time 


It is of interest also to consider the spacewi.se description of the 


process. In this case the time lag must be replaced by the 


le Lai 


which is a vectorial quantity indicating the location in the duct where 

the transformation from the "heavy" to the "light" fluid takes place. 

The space lag is denoted by ^ on Figure 11=3. Of course 3 the space 

lag can be related to the time lag when the particle velocity is known. 

Like in combustions the location in the duct where the transformation 

„ 

takes place can be regarded as the source of the light fluid. It is 

» it 

obvious that the flow properties in the region occupied by the light 
fluid will depend upon the intensity of this source. If it is assumed 

ii, ii 

that the injection rate of the heavy fluid is constant and that the time 
and space lags were constant then the intensity of the source would also 
be independent of time resulting in a .steady flow in the "light" fluid 
region. However 3 this is not the case because fluctuations which affect 


the time lag and/or the injection rate are present both at supercritical 
and at subcritical pressures. In the vicinity of the critical thermo- 
dynamic point large fluctuations of properties , in particular of density, 
are observed even in non- flow systems . In boiling systems fluctuations 
are always present because of variations of the rate of bubble formation 
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! FIGURE 1-3 "TWO-REGION" APPROXIMATION SHOWING THE 

TIME LAG ANO THE SPACE LAG. 


and population, of flow regimes, of the heat transfer coefficient, etc. 
Consequently, the strength of the source may fluctuate even when the in- 
jection rate is kept constant. It is evident also that variations of 
inlet velocity will introduce additional effects. 

The nucleation and evaporation at subcritical pressures and the 
transformation of "heavy" clusters to "light" clusters at supercritical 
pressures are rate processes that occur during and have an effect upon 
the length of the time lag. Both of these transformation rates are af- 
fected by the pressure, temperature and by other rate processes such as 
the rate of energy transfer, flow rate, et.c. If one of these factors 
changes or fluctuates, the transformation rates will fluctuate also 
resulting in a fluctuation of the time lag, i.e», in the fluctuation 
of the source. Since the source affects the flow conditions in the "light" 
fluid region the flow in this region may become oscillatory. 

II. 3 Organized Oscillations 

Oscillations of a system can be always produced if properly excited. 

Such oscillations can be distinguished by a characteristic time,. i. 6-., 
period if the process is periodic or by a relaxation time if it is 
aperiodic. 

Like in combustion and following Crocco and Cheng M we shall 
distinguish two cases: random fluctuations and organized or coordinated 

oscillations. 

As random fluctuation we consider those that are Similar, to fluctuations 
observed in ordinary turbulent flow. In this case it can be assumed that 
the transformation process, for example the rate of evaporation in boiling 
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at subcritieal pressures, is not excited. The fluctuations at one point 
do not have any effect on other fluctuations somewhere else in the system. 

n 

Since the integrated effects of these fluctuations vanish they do not pose 
a problem. 

In the case of an organized oscillation the transformation process 
will be excited by one or more coordinating processes such as the oscillation 
of the inlet flow rate, of the heat transfer coefficients etc. The exciting 
force for maintaining the oscillation of the coordinating process is in 
turn provided by the transformation process. For examples in boiling sys- 
tems oscillations of pressure will affect the saturation temperature which 
may induce oscillations in the rates of evaporation. These in turn may 
induce flow oscillations and provide the excitation force for maintaining 
the pressure fluctuations. 

The fundamental character of organized oscillations is that a well 
defined correlation exists between fluctuations at two different points 
or instants. In other words that a disturbance is propagated., i.e. s dis- 
placed in time and space through the system. When these organized oscil- 
lations are present their integrated effect does not vanish whence the interest 
in these oscillations. Furthermore , an oscillatory system may become unstable „ 
i.e. s it may have the tendency to amplify. In the example cited above 
pressure fluctuations of an increasing amplitude may be generated leading 
to the destruction of the system. When the effects are proportional to 
the causes the system is defined as linear. In this paper we are interested 
in such systems. 
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I I . 4 The Mechanism of Low Frequency (Chugging) Oscillation at 
Supercritical Pressures 


It was noted in the preceding section that the characteristic of 
organized oscillations is the propagation of disturbances through one system. 
These disturbances can be variations of density, pressure, enthalpy, entropy, 
etc. In this section we shall examine the effects which these propagations 
may have on the oscillating propensity of the system. In particular, we 
shall consider the propagation of density disturbances and the effect of the 
time lag, i.e. , of the space lag. The effects of pressure waves are discussed 
in Section 11-7 together with the other mechanisms which may induce flow 
oscillations in the nearcritical and supercritical regions. 

We note that the effect of the time lag in inducing combustion in- 
stabilities was already analyzed by Summerfeld [52^ , Crocco and Cheng [ 48^ 
among others. In boiling systems, this effect was already analyzed by 
Frofos [ 49 , Wallis and Heasley 5o""j and by Bour& [4 In these analyses 

the flow was assumed to be homogeneous, i.e., the effect of the relative 
velocity between the gas the liquid phase was neglected. ,A density propa- 
gation equation, applicable to two-phase mixtures, which takes this effect 
into account was formulated in [53~| and solved in [.54, 55^ . 

Let:, us examine now the effects of the finite rates of propagation and 
of the resulting time lag and time delays on the flow in a system consisting 
of a constant pressure tank connected by a feeding system to the heated duct. 

Consider first tjae tank and the feeding system only and let us perturb 
suddenly the inlet flow. If there is no feedback between the heated unit 
and the y.pstream part of the system, the steady state conditions will be 
restored. In particular, if the variation of the flow rate is small during 
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FIGURE 1-4 THE EFFECTS OF TIME LAG AND OF SPACE LAG 
IN INDUCING FLOW OSCILLATION. 








the time a pressure wave propagates back and forth through the tank and 
the feed system, then the pressure effects can be neglected. As discussed 
in £ 48^ the process can be described then, with sufficient accuracy, by 
an exponentially decaying flow which is characterized by the relaxation 
time constant i.e., by the line relaxation time. Therefore, the system 
is stable because the steady state conditions will be eventually restored,* 

Consider now the effect of a perturbed flow at the inlet of the heated 
duct (See Figure II-4), It is obvious that an oscillatory flow at the inlet 
will induce an oscillatory flow of the fluid in the duct. However, in 
absence of a driving mechanism these oscillations would also exhibit an ex- 
ponential decay. We are looking for a mechanism whereby these flow oscil- 
lations at supercritical pressures can be maintained. Like in boiling and 
in combustion such a mechanism is provided by the propagation phenomena 
which introduce different delay times in the response of the system. This 
is shown in Figure II-4. 

It can be seen on Figure II-4 that an oscillatory inlet flow can induce 
oscillations of the space lag; this is in accordance with the discussion of 
the preceding section. The onset of these oscillations is delayed however 
by the lag time because of the finite rate of propagation of the dis- 

turbance, An oscillatory space lag, which is equivalent to an oscillatory 
source, will induce flow oscillation, in the "light" fluid region. These 
source- induced oscillations will be present in addition to those already 
induced by the inlet flow. Because of these two oscillatory motions there 
will be a delay time ©•* in the flow response. Oscillations of. the flow 

*We are assuming here that any servo-control mechanism in the feeding system 
will not have a destabilizing effect, 
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will induce oscillations of enthalpy and of density both delayed by a 
certain delay time. With flow and density oscillating, the pressure drop 
in the duct will also oscillate. If the conditions are such that the mini- 


mum pressure drop in the duct occurs when the inlet flow is maximum, it is 
apparent then that the oscillations can be maintained. It is also obvious 
that whether or not this will occur will depend on the time lag ~C ^ and on 
the delay times . When these delay times do not 

depend upon H it can be seen that increasing the lag time X ^ has a 
destabilizing effect. Since the time lag X*, (see Figure II- 4) depends 
upon the enthalpy difference i - i^, it can be concluded that, for this 
particular case, a decrease of inlet enthalpy i^, i.e. , that an increase 
of ^^21 k aS a destabilizing effect. 

From this qualitative description it can be already seen that at 
supercritical pressures an unstable flow can be induced by the delayed 
response of various perturbations. It remains now to advance a qualitative 
description. We shall do this in the following chapters by modifying and 
applying the method proposed in ^50, 5l”^ for boiling at subcritical 


pressures . 

II . 5 Method of Solution 

In what follows we shall consider the "heavy" and the "light" fluid 
regions separately. Each will be described in terms of three conservation 
equations and the equation of state. We shall use the one dimensional 
representation and obtain solutions for each region. These solutions will 
be matched at the transition point, i.e., at the end point of the space 
lag to provide a solution valid for the entire system. 


Following ^48. 50, 5lJ it will be assumed that the variation of 
pressures can be neglected,, This is implied by the assumption that the 
density is function of enthalpy only . It can be seen that this assumption 
will be valid only if the variations of flow, density, enthalpy, etc. are 
relatively small during the total time for propagation of a pressure wave 
back and forth through the duct. Under this condition it can be assumed 
that the various disturbances move through a uniform medium. It is ap- 
parent also that this will be true only if the rate of propagation of 
pressure waves is considerably faster than the rate of: propagation of the 
disturbances. However, both in boiling systems as well as in the nearcritical 
region the velocity of sound reaches ve t; ' low values.* Consequently, it can 
be expected that there will be a range of operating conditions for which 
the assumption that the properties do not depend upon pressure variations 
will not be satisfied. For boiling systems this limitation has been already 
recognized and discussed by Christensen and Solberg ^56~j „ In general, 

it can be expected that the assumption will be satisfied in the low frequency 
range, i.e., in '‘chugging” oscillations. When the effects of pressure vari- 
ations can be neglected then one can use the formulation put forward in £50~) 
and carried out in (j>l^ for boiling systems at subcritical pressures. 

The method of solution used in this paper is as follows. A small per- 
turbation is imparted to the inlet velocity. The velocity of the fluid is 
determined then by integrating the divergence of the velocity. With the 


^Indeed, in the critical region some authors reported values approaching a 
zero velocity. At present neither the exact value of the velocity of sound 
at the critical point is available nor a satisfactory understanding of the 
phenomenon has been attained* 
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velocity known the energy equation is integrated to obtain the time lag X. ^ 
as well as the rate of propagation of enthalpy disturbances., From the 
enthalpy and from the equation of state we then obtain the density of the 
medium. The differentiation between the nearcritical region and the 
supercritical region is achieved by assigning the appropriate expression 
to the equation of state. With the velocity and the density known, the 
momentum equation can be integrated. Since the inlet disturbance is small, 
the momentum equation is first linearized and then integrated to give the 
characteristic equation. 

Because of the linearization of the momentum equation the analysis 
will be applicable only to cases where the effects of the instability are 
not so strong to produce large amplitude oscillations. It can be used there- 
fore to predict the conditions of incipient instability, i.e., to determine 
stability limits. As discussed in [48 ] and ^50, 51^ linear effects and formu- 
lations have been successful in predicting certain type of instabilities 
("chugging" instabilities) in combustion and in boiling systems respectively. 

A similar result could be expected, therefore, with the present formulation if 
it is used to predict the onset of "chugging" instabilities at supercritical 
pressures . 

■*a> 

II. 6 The Characteristic Equation and Its Applications 

The characteristic equation for this problem is an exponential poly- 
nomial, it is therefore of the same form as the characteristic equation for 
combustion instabilities ^48 \ , thus 

(U) , L,d) — e ST * t a (s) = o n-i 
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where and L ^ are polynomials with coefficients independent of the time 
and where s is a root of the characteristic equation. 

In general, the root s is a complex number; the real part gives the 
amplification coefficient of the particular oscillatory mode, whereas the 
imaginary part represents the angular frequency. Since the original per- 
turbation is assumed to be of the form expj^ st , a given oscillatory mode 
will be stable, neutral or unstable depending upon whether the real part 
of s is less, equal or greater than zero. A sufficient condition for the 
system to be stable is therefore that the characteristic equation (Eq. II-l) 
has no roots in the right half of the complex s plane. 

Let us examine now what information can be obtained from the character- 
istic equation as well as the type of practical problems where this information 
can be applied most usefully. Two such problems were discussed by drocco and 
Cheng ^ 48"] in connection with the stability analysis of combustion systems. 

The same discussion can be applied to the present problem. 

In the first class of practical problems one is interested irideter- 
mining whether a given system with specified characteristics, i.e., with 
specified numerical coefficients is stable or unstable. This is most often 
a situation that arises during the planning period; i.e., before the system 
is designed and tested. The characteristic equation can be used to provide 
an answer to this type of problem. In particular, since the numerical co- 
efficients in the characteristic equation are known, Crocco and Cheng £48<”j 
note that the use of Satche's [58^ diagram is most useful for analyzing the 
exponential polinomial obtaining thereby a solution for this type of problem. 

In the second class of practical problems one is interested in the 
qualitative trends of the stability behaviour of a system when various 
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parameters are changed. This is most often a situation that arises during 
the design period because of the designer’s need either to design a system 
with sufficient safety margins or to modify a given unstable system in order 
to make it stable. For this kind of problem Crocco and Cheng £48~~j 
note that it is advantageous to use the characteristic equation for deter- 
mining the stab! ity boundary of a certain system. On such a boundary, 
expressed in terms of the operating characteristics of the system and of 
the process, the oscillatory mode in question is neither stable nor unstable, 
i.e., the real part of s vanishes for that mode. The stability boundary 
divides therefore the space formed by the parameters of a given system into 
different domains in which the system is stable on one side of the boundary 
and unstable on the other. If by varying one parameter of the system the 
stability boundary is shifted in such a way as to decrease the unstable 
domain, the variation of the parameter has a stabilizing effect and vice 
versa. 

Following the standard procedure the stability boundary is obtained 
from the characteristic equation by setting s = i &-> , where to is the 
frequency of the neutral oscillation. Upon separating the real and imaginary 
parts of the characteristic equation one can eliminate the frequency Co , 
the resulting equation represents the stability boundary. Two such boundaries, 
obtained from characteristic equations given by first and second order ex- 
ponential polynomials, are shown in subsequent chapters. 

II. 7 Additional Mechanisms Leading to Unstable Operation 

Before proceeding with the formulation of the present problem we shall 
note and examine briefly additional mechanisms which can induce flow oscillations 
in the nearcritical and supercritical regions. These mechanisms will be 
analyzed in more detail in separate publications. 


- 25 - 



It is instructive to note first a general characteristic of oscillatory 
systems. A necessary condition for maintaining oscillations is that enough 
energy is supplied to the system at the proper frequency and phase relation 
in order to overcome the losses due to various damping effects which are 
always present in real systems. When the rate of energy supplied is control- 
led by an external source and is independent of the fluctuations inside the 
systems, the oscillations will build up when the energy is released at a 
characteristic frequency giving rise to the resonance phenomenon. However, 
when the system contains itself an energy source, with a property that the 
energy release depends upon a fluctuation inside the system, then an accidental 
small disturbance in the system may interact with the source resulting in 
oscillations of increasing amplitude. For this to take place it is necessary 
that the energy from the source be fed to the disturbance during part of 
the cycle. 

It was discussed in preceding sections that the system which is analyzed 
in this paper has the property that the energy release depends upon fluctu- 
ations inside the system. Oscillations of the time lag and of the space lag 
are examples of such fluctuations. We shall examine now other energy sources 
and fluctuations which may be present in the system. 

It 

oscillations can occur in a system consisting of a gas flowing through a 
heated duct. For such oscillations to be maintained Rayleigh notes that 
the energy must be added to the gas during the moment of greatest conden- 
sation and removed during the rarefaction period. This leads to the Rayleigh's 
criterion which states that a component of the fluctuating heat addition must 
bq in phase with the pressure wave if oscillations are to be thermally driven. 
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was discussed by Rayleigh |~59”| that internally driven pressure 


The same criterion can be used to explain a type of oscillation ob- 
served at critical and supercritical pressures as well as in boiling 
mixtures at subcritical pressures. In both systems the heat transfer co- 
efficient is a strong function of pressure. Thus, pressure oscillations 
may interact with the heat transfer coefficient inducing oscillations of 
the latter. If these oscillations are in phase, the system may be thermally 
driven and become oscillatory. 

Another mechanism which may induce oscillations at both subcritical 
and at supercritical pressures is caused by the large compressibility of 
some parts of the system. At high pressures this is the section of the 
system where the properties of the fluid pass through the nearcritical 
region. At pressures below the critical point, this will be the section 
of the system where subcooled boiling takes place. 

Still another mechanism that can induce oscillations at subcritical 
pressures is caused by the change of flow regime which can induce large 
fluctuation of the mixture density. These in turn may induce both os- 
cillatipns of the flow and of the heat transfer coefficient thus providing 
the driving force necessary for maintaining the oscillations. 

It can be expected that each of the mechanises may be effective over 

i 

some range of operating conditions. It can be also expected that the 
resulting oscillations are characterized by a certain frequency range and 
by particular wave forms. Indeed, several frequency ranges and wave shapes 
have been reported and described in the references discussed in Chapter I. 

The mechanisms just described will be the subject of future investigations; 
in what follows we shall proceed with an analyses of the "chugging" oscillations. 
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III. The “Heavy" Fluid Region 


III.l The Governing Equation 

For a "two-region" approximation the "heavy" fluid region extends 
from the entrance of the heated duct to the transition point. Note, 
that for a system with constant energy input, the location of this point 
will move along the duct when the inlet velocity and/or the inlet enthalpy 
vary. 

It will be assumed that the fluid in this region is incompressible 
and that the thermodynamic properties are constant. The problem is formu- 
lated by considering the three field equations describing the conservation 
of mass, momentum and of energy in addition to the constitutive equation 
of state describing the properties of the fluid. 

For a one dimensional formulation, used in this analysis, the con- 
tinuity, energy and momentum equations are given respectively by: 
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where the symbols are defined in the. Nomenclature. 
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The constitutive equation of state is given by: 



C°* it, 


III-4 


Equations III-l, 2, 3 and 4 are four equations which specify the four 

variables ^ and i in the "heavy" fluid region. These four 

equations will be integrated to yield T. P 'M and i as function of space 

' ? 

and of time. These will be then used to determine the time lag and the 
space lag. 

III. 2 The Equation of Continuity and the Divergence of the Velocity 

In view of the assumption of an incompressible "heavy" fluid the con- 
tinuity equation reduces to the divergence of the velocity 



^ HI “5 

whence upon integration we obtain 

Us U (t) III-6 

The velocity in the "heavy" fluid region is therefore independent of 
position, it is a function of time only. 

In order to analyze the stability problem we shall assume that a small, 
time dependent velocity variation Sm s , is superimposed on a steady inlet 
velocity , thus: 


where the bar indicates steady state conditions., 


III-3 The Energy Equation 

With the fluid velocity given by Eq. III-7, the energy equation becomes 
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This is a first order partial differential equation whose solution can be 
obtained by merxs of characteristics ^60, 6l] . The general solution of 
Eq. III-8 is of the form; 


t = f (t) 
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where 
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are solutions of any two independent differential equations which imply 
the relationships: 
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For example, by taking alternately the first and the second equation, the 
first and the third equation we obtain the following set 
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In order to solve the problem it is necessary to specify the initial 
and the boundary conditions. These will be specified by letting a particle 
enter the duct with enthalpy i^ at time X, , (See Figure II-3) thus, 
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With this boundary and initial condition, one obtains after integrating 
Eq. III-12 and III-13 the following relations: 
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Upon eliminating the time t“T ( > between these two equations we obtain 
an expression for the enthalpy as function of space and time, thus: 
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The first term on the right-hand side is the steady state term, whereas 


the second one is the transient which accounts for the perturbation of the 
inlet velocity. 


III. 4 The Time Lag and The Space Lag 

In Section II- 2 the time lag TT^, was defined as the time required for 
increasing the enthalpy of the fluid from the inlet value i^ up to the 
enthalpy at the transition point i (See Figure II-3). Consider now a 
fluid which enters the duct at time TTj and attains the enthalpy i^ at 
time ~C 2 » it follows then from Eq. III-X6 that the time lag is given by: 
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which, in view of Eq. Ill- 13 can be also expressed as: 
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It can be seen that for a given system and at a given pressure the time 
lag depends only upon the enthalpy difference (i ^ “ i^ * A-i-t\ ) an ^ the 
heat flux density. 

We shall determine now the space lag. For a "two region" approximation. 
Figure II-3 indicates that the "heavy "fluid region extends up to the transition 
point where the bulk fluid enthalpy reaches the value of i 0 . Inserting i ^ 
in Eq. III-17 we obtain the following expression for the space lag: 
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For steady state operation £. = 0, whence we obtain from Eq. III-2Q 
the steady state space lag /\ , thus: 
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In view of Eq. III-18 and III-21, we can also express Eq. III-20 as: 

St “ 

Aft) = -t- 

or 

A (f) =5 A + $ A Hi-23 

Several comments are appropriate. 

1) Equations III-18 and III-21 indicate that for steady state, i.e., 

when £, = 0 the time lag corresponds to the transit time. of a fluid 

particle through the "heavy" fluid region. 

2) Equation III-22 shows that the response of the space lag to a 
variation of the inlet velocity is delayed by a time period equal to the 
time lag ”£ b . 

3) If we interpret the enthalpy i ^ as the enthalpy at saturation and 

therefore the difference ^ i the subcooling then Eq. III-22 predicts 

the location of the boiling boundary, i.e., the location where boiling 
starts in a two-phase mixture at subcritical pressures. Indeed, such an 
expression was derived previously (in [49, 50, 31 J among others, using 
somewhat different approaches) for analyzipg oscillations in boiling mixtures 
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The time lag X^ was called there the "evaporation time constant" [49j . 
Ill .5 The Momentum Equation 

The momentum equation can be integrated now since the velocity in 
and the boundaries of the "heavy" fluid region have been determined. With 
the boundary conditions taken as 


T* T>, 
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the integrated momentum equation becomes 
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where we have taken into account the assumption that the density in the 
ri heavy H f luid region is constant. In view of Eq. III-5, III-6, III-7, and 
III-23, the integrated momentum equation yields: 
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Linearizing Eq. III-26 and retaining only the terms with the first power 
of £. we obtain: 
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We shall consider now the pressure drop across the inlet orifice. 
Defining by a numerical coefficient that takes into account the effect 
of the geometry of the restriction and of other losses like vena contracta, 
etc. , we can express the inlet pressure drop across the inlet orifice by: 
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which, upon linearization can be expressed as: 
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We define now the steady state values of the pressure drop due to 
body forces (gravity) by: 
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due to friction by: 
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and due to the inlet orifice by; 
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In view of these three relations and upon substituting Eq. III-29 
in Eq. III-27, we can express the pressure drop in the "heavy" fluid 
region by; 

— ?o = A?o, + 


r f a dsu, 

^ ?01 

+ & M. 



3a?,J 

H 1 at 

1 


1 

ll 

OK J 


III-33 


where 


$ M, 


£ e 


st 


III-34 


and 


; S~t 

<TA = (i 

S 


-SZl 

e ) 


III-35 


The first line on tile right-hand side of Eq» v III-33 represents 
the sura of the steady state pressure drops, whereas the second one accounts 
for the transient response. In particular, the first terra is the inertia 
of the "heavy" fluid region; the second termaire the pressure losses due to 
variation of inlet velocity whereas the lhst term shows the effect of the 
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varying space lag. Equation III-35 indicates that this last effect is 
delayed by time lag 71^. We shall proceed now with the analysis of the 
"light" fluid region. 



IV The "Light 11 Fluid Region 


IV . 1 The Governing Equation 

For a "two region" approximation the "light" fluid region extends 
from the transition point to the exit of the heated duct. The problem 
is formulated again in terms of three field conservation equations and 
of a constitutive equation of state. However, in contrast to the "heavy" 
fluid the density in the "light" fluid region is function of enthalpy 
and of pressure. It was discussed in Chapter II that for "chugging" 
oscillations, the effects of pressure variations can be neglected. 
Consequently, the density will be a function of enthalpy only. 

The "light" fluid region is described, therefore, in terms of the 
continuity equation 




'Dt 




2i 

% 




Oa 




n 


= o 


IV -1 


of the energy equation 
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the momentum equation 
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and the constitutive equation of state 
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or 

or — vii) iv-5 

when expressed in terms of the specific volume ’’0* . 

Equations IV- 1, 2, 3 and 4 are four equations in terms of four 
variables *P, ^ , 'U. and i. They are applicable to the "light" fluid 
region at supercritical pressures. They can be also applied to the 
two phase flow region at subcritical pressures if, and only if, the 
relative velocity between the phases can be neglected. 

It is emphasized here that the form of the energy and the form of 
the momentum equation, as given by Eq. IV"2 and Eq. IV-3, are correct 
only if the relative velocity between the phases is either zero or its 
effect is negligibly small. If such is not the case, then both Eq. IV-2, 
and Eq. IV-3 must be modified . 

It was discussed in Section II-5 that at high pressures, say above 
0.85 of the reduced pressures, the effect of the relative velocity is 
so small that it can be neglected. The region of interest to this analysis 
is the high pressure region. It can be expected, therefore, that, in 
this investigation, both Eq. IV-2 and Eq. IV-3 can be used to approximate, 
with sufficient accuracy, the energy and the momentum equation for the 
two phase mixture in the nearcritical region. 

In what follows we shall use, therefore, Eq. IV- 1 through 4 to des- 
cribe both the "light" fluid at supercritical pressures and the two phase 
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mixture in the nearcritical region. The differentiation between the 
"light" fluid and the two phase mixture will be realized by assigning 
the appropriate expression to the constitutive equation of state. This 
will be done in the section that follows. 


IV. 2 The Equation of State 

For the "light" fluid the relation between the specific volume and 
the enthalpy can be obtained either empirically, i.e., from experimental 
data or it can be approximated by an equation of state such as the per- 
feet gas or the van der Waals ' equation etc. It was noted in Section II. 1 
that for oxygen the perfect gas equation predicts with sufficient accuracy 
the relation between the specific volume and the enthalpy. Since this 
fluid is of primary interest to this investigation, the perfect gas 
equation will be used as the constitutive equation of state for the "light" 
fluid at supercritical pressures. 

Assuming a constant pressure process we have for a perfect gas the 
following relations 
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For a "two-region" approximation the boundary condition for the "light' 
fluid region is given by: 

V - Vj. at i , i t IV _ 9 

* 

We obtain then the equation of state for the "light" fluid region by 
integrating Eq. IV-8 subject to the boundary condition given by Eq. IV-8, 
thus: 
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We shall derive now the equation of state for the two phase mix- 
ture in the nearcritical region. We recall first that the quality x, 
of a two phase mixture is defined by: 
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where G and G f are the mass flow rates of the vapor phase and of the 
§ >t 

liquid phase respectively. We recall also that the specific volume 
and the enthalpy of a two-phase mixture are given by: 


v - (i-X)Vj. +- A vr 
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Where V - „ i and IT , i are the specific volume and the enthalpy of 
f f g g 

the liquid and of the vapor respectively. 

We obtain then the equation of state for the two phase mixture by 
eliminating the quality x, between Eq. IV-12 and Eq. IV-13, thus: 
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Where ^ 'V = V 1 - V , and where ^ i is the latent heat of vapori- 
tg g t tg 

zation. Differentiating Eq. IV- 14 we obtain for the two phase mixture: 



IV- 15 


which can be compared to Eq. IV- 8 applicable to the "light" fluid at 
supercritical pressures. 

It is important to note that both, the equations of state for the 
"light" fluid at supercritical pressures., i.e. , Eq. IV-10, and the equation 
of state for a two-phase mixture at subcritical pressures., i„e., Eq. IV-14, 
are of the same form, i.e., both can be expressed as: 

oJCt) = ^ -v- iv- 16 

We can use, therefore, Eq. IV-16 for the equation of state in the near- 
critical as well as in the supercritical region. We shall distinguish 
one region from the other by substituting either Eq. IV-15 or Eq. IV-8 
in Eq. IV-16. 
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IV ,3 The Equation of Continuity and the Divergence of Velocit' 


Several methods are available ^49, 50, 51, 62^ for determining the 
velocity in a boiling mixture. Any of these could be modified and used 
to determine the velocity in the "light" fluid region. In what follows 


we shall use the method of Bour£ 
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As in the "heavy" fluid region we shall determine the velocity by 
integrating the divergence. However, in contrast to the "heavy" fluid 
region where the divergence is given by Eq. III-5, the divergence in the 
"light" fluid region is not zero but is obtained from Eq 0 IV- 1, thus: 
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In order to integrate the divergence it is necessary to evaluate 
the right-hand side of Eq. IV-17. Following Bourl this cAn be done by 
means of the energy equation. 

Since the density is function of enthalpy only one can write 
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Substituting Eq. IV- 2 in Eq. IV-18, it follows that: 
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whence from Eq. IV-19 and IV-17 we can express the divergence as: 
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We shall define now the reaction frequency* SI by: ' 
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It follows then from Eq. IV-8 and Eq. IV-21 that the reaction frequency 
for the "light’' fluid in the supercritical region is given by: 
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whereas it follows from Eq. IV-15 and Eq. IV-21 that for boiling at sub 
critical pressures the reaction frequency is given by: . 
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With the reaction frequency defined by Eq. IV-21, it follows then 
from Eq. IV-20 that we can express the divergence as: 
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*The reasons for using this definition are discussed in Section IV- 9 


The physical significance of this equation is simple” the divergence 
of the velocity in the "light" fluid region is equal to the volumetric 
rate of formation of the "light" fluid per unit volume of space. 

In order to integrate Eq. IV-24 it is necessary to specify the 
boundary conditions; these are given by considering the velocity in the 
"heavy" fluid region 3 i.e. 9 Eq. III-7. The boundary condition for 
Eq. IV-24 is therefore? 




St 


at 7 w X It) 

IV- 25 


whence upon integration of Eq, IV-24 , we obtain for the velocity in the 
"light" fluid region the following expression? 

uiiitj = U, t e.e A(t)"] IV-26 

We note that Eq, IV-26 with _TL given by Eq. IV-23 is the velocity 

in the two phase boiling mixtures as such it was used already in (49,50, 51and 62.) 

It is instructive to examine further Eq. IV-26 S which 9 * by means of 
Eq. III-23s can be expressed as? 

it ^ ,-S C 
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jr 

We obtain the, steady state velocity in the "light" fluid region by letting 
= 0 in Eq. IV-27* thus? 
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We can rewrite then Eq. IV- 2 7 as: 

Ul\,t) » Uyi^j + f^lt) IV- 2 9 

where the time dependent perturbation of the "light" fluid is given by: 
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( It can be seen from Eq. IV-29 and IV-30 that the flow in the "light 1 
fluid region is affected by both the inlet perturbation as well as the 
perturbation of the space lag. This last perturbation is delayed by the 
time lag X. (see Figure II-4). 

We shall define now several steady state relations which shall be 
used in following chapters. 

By letting z = -£■ in Eq. IV-28 we obtain the steady state velocity 
at the exit of the heated duct, thus: 


The lengthwise average velocity in the "light" fluid region is defined 
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whence from Eq. IV- 28 we obtain: 
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From Eq. IV- 31 we have: 
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Substituting this relation in Eq. IV-33 we have the following expressions 
for the average velocity: 
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The log mean velocity in the "light" fluid region is defined by: 
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where we have taken into account Eq. IV- 34. 
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A fourth relation of interest to this investigation can be obtained 
from the conservation of momentum G*and the definition of the log mean 
density,, If we denote by ^ the density of the fluid at the exit from 
the heated duct, then the log mean density in the "light" fluid region is 
given by: 
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The mean velocity , based on the log mean density, is then obtained by 
considering the mass flux density, i.e., the momentum G, thus: 

G 

r~ = e IV 

which, in view of the preceding relations can be expressed also as: 
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We can proceed now with the solution or the energy equation. 

IV-4 The Energy Equation 

We can solve the Energy equation now since the Velocity and the equation 
of state in the "light" fluid region are specified. By substituting Eq. IV- 26 
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and IV-16 in Eq. IV-2 we can express the energy equation as: 
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Taking the enthalpy i ^ at the transition point for reference and in view 
of the definition of the reaction frequency JH. , given by Eq. IV-21, 
we can rewrite Eq. IV-36 as: 
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The initial and boundary conditions for the energy equation are de- 
termined by considering the conditions at the transition point (See Figure 
II-3); they are given therefore by: 

Ul, = 0 crf-.t.-Ti. 
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i-k = 0 ■-* A ♦ £.e 

IV- 43 

Equation IV-41 is again a linear partial differential equation, it 
can be solved therefore by the method used in solving the energy equation 
for the "heavy" fluid (See Section III-3). Following this procedure, we 
obtain in place of Eq. III-ll the following set: 
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Integrating first Eq, IV-46 and taking into account the initial con' 
ditions given by Eq. IV-42 we obtain: 
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whence : 




In order to integrate Eq. IV-45 we note that it is of the form of: 
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whose solution is: 
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The integral of Eq. IV-45 is therefore: 
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which, in view of Eq. IV-48, can be expressed as: 
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The valtie of the constant of integration is evaluated by means of Eq. IV-42 
and IV-43, thus: 
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Substituting Eq. IV-53 in Eq. IV-51 we obtain after some re- 


arrangement: 
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which, in view of Eq. III-23 and Eq. IV-30, can be expressed as 
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By substituting Eq. IV-48 in Eq. IV-55 we obtain the solution of the 
energy equation for the specified initial and boundary conditions , thus 
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from which we obtain the enthalpy for the "light” fluid region 
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Expanding and retaining only the first power of 6 we obtain after some 
rearrangement 
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If we let the perturbation go to zero, i.e., £ = 0, we obtain from 

Eq. IV-58 the enthalpy for steady state operation, thus 
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IV- 5 The Residence Time 

It is of interest to evaluate now the steady state residence, i.e., 
the transit time of a particle in the heated duct. Denoting by ig the 
enthalpy at the exit and by 'C 3 the time when a particle reaches this 
enthalpy we obtain from Eq. IV-47 the residence time in the "light" 
fluid region, thus : 
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which, in view of Eq. IV-21, can be expressed also as: 
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Denoting by Q = the rate of energy transfer to the entire 

duct and by , the mass flow rate, we can express the total energy 


balance in steady state as: 
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Substituting this relation in Eq. IV-61 and in view of Eq. III-18 we 
obtain the following expression for the residence time in the heated 
duct : 
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IV-6 The Density and the Density Perturbation 

The density in the "light" fluid region is given by the equation of 
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state, i.e., by Eq. IV-16, which, in view of the definition of the re- 
action frequency _TL , given by Eq. IV-21, can be expressed also as: 
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Since the enthalpy in the "light" fluid region is given by the 
solution of the energy equation, i.e., by Eq. IV-56 we can express the 
density as function of time or as function of time and space. Thus by 
substituting Eq. IV-64- in Eq. IV-47 we obtain: 
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whereas by substituting Eq. IV-64 in Eq, IV-56 we obtain: 
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which, by means of Eq. IV-30, III-20 and III-18, can be expressed also 


or, in view of Eq. IV-65, it can be transformed in 
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Using again the binomial expansion and retaining only the first power 


of £ 

we can 

express 

Eq. IV- 

•67 as: 
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Similarly, we 

can expand Eq. 
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where the steady state velocity U^(^) is given by Eq. IV-28. 

By letting ^ = 0 in Eq. IV-66 we can obtain the local steady state 
density in the "light" fluid region, thus: 
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We shall define now several steady state relations which will be 


used in the following sections. 

By letting ^ in Eq. IV-71 we obtain the density 

from the heated duct thus : 
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where we have taken into account Eq. IV-31. 

We shall define now the average density in the "light" 

by: 
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whence from Eq. IV-71 we obtain: 
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In view of the definition of the log mean velocity give 
this average density can be expressed as: 
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by Eq. IV-36 
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We have already defined the log mean density by: 
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where ^ i s given by Eq. IV-72. 

A fourth expression can be obtained from the definition of the 
average velocity given by Eq. IV-35 and the momentum (3r . We 
can express therefore a mean density , based on the average velocity, 
by: 
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With the steady state density in the "light" fluid vapor given by 
Eq. I.V-71, we can express Eq. IV-69 as 
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where the density perturbation is given by 



which, in view of Eq. IV-30, can be expressed also as: 
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By letting ^ ^ in Eq. IV-79 we obtain the density perturbation 
at the exit from the heated duct, thus 
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It can be seen from the preceeding equations that in the "light" 
fluid region the density perturbation is affected by both the perturbation 
of the inlet velocity and by the variation of the space lag. Further- 
more, the effect of the inlet velocity perturbation is delayed by a 
delay time. Equations IV-30 and IV-80 are the quantitative expressions 
for the flow and density variations in the "light" fluid region whiph 
were qualitatively described in Section II r4. 
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With the density and the velocity in the "light" fluid region 
given by the expressions derived in this section and in Section IV-3 
respectively, we are in the position to integrate the momentum equation. 
IV. 7 The Momentum Equation 

In order to integrate the momentum equation it is necessary to 
specify the boundary conditions, these are given by: 


?= \ 





or f 

i - 1 

IV- 82 


whence the integrated momentum equation becomes: 
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The expressions for the density and the velocity which should be 
substituted in this equation are given by Eq. IV-78 and Eq. IV-27, 
i.e., Eq. IV-29 respectively. We shall consider now each term of Eq. 
IV-83 separately. 

IV .7 . 1 The Inertia Term 

The inertia term in the momentum equation is given by: 
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Substituting Eq„ IV-78 and Eq. IV-27 in Eq. IV-84 and retaining only 
the first power in £. we get: 
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In view of the definition of the average density and of the velocity 
perturbation given by Eq. IV-75 and Eq. IV-30, respectively, the inertia 
term can be expressed as: 
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IV . 7 . 2 The Convective Acceleration Term 

The convective acceleration term in Eq. IV-83 is given by: 
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Substituting Eq. IV-78 and Eq. IV-27 in Eq. IV-87 and retaining only the 
first power in ^ we obtain upon integration: 
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It is of interest to examine the physical significance of the various 
terms . 

If we let * 2 , = 0 in Eq. IV-88, we obtain the steady state acceleration 
pressure drop 4P« , which, in view of the definitions given in Sections 

IV-3 and IV-6, can be expressed as: 

A?« = 6 nU-A) .^'(cr-j-u,) = <f,'jlt^(Mj-W 1 ) = 

Ui 

— .-U*, ( M, ) 

The second term in Eq. IV-88 can be expressed by means of the 
Eq. IV- 89 and of the space lag variation defined by Eq. III-22, thus 
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It shows, therefore, the influence of the variation of the space lag on 
the acceleration pressure drop in the "light" fluid region. 

In view of Eq. IV-89 and Eq. IV-30, the third term in Eq. IV-88 
can be expressed as 
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It expresses, therefore, the influence of the velocity perturbation in 
the "light" fluid region on the acceleration pressure drop. 

The last two terms in Eq. IV-88 stem from the density perturbation 
term in Eq. IV-70, i.e., from 
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The third term in Eq. IV-88, i.e., the first term on the right hand side 
of Eq. IV-92 can be expressed in terms of Eq. IV-89 thus 
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It shows , therefore, the effect, of the variation of the velocity in the 
"light" fluid region on the density and, therefore, on the acceleration 
pressure drop in that region. 

The physical meaning of the last term on Eq. IV-88, i.e., Eq. IV-92 
is not as clear as that of the other terms in Eq. IV-88. An insight can 
be gained however, by considering the upper and lower limits of the 
integral 
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It is shown in the Appendix C that this integral is bounded by: 
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which, in view of Eq. IV-89, can be expressed as 
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We note that a 
in Eq. IV-94, 
the integral I 


simple expression can be obtained by setting 
this approximation results in the following expression for 
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The physical meaning of the integral 1^ is now clear: it expresses the 
effect of the perturbation of the inlet velocity on the density (see 
Eq. IV-69) and, therefore, on the acceleration pressure drop in the 
"light" fluid region. This effect is delayed by a delay time equal to 
"Z\,or to ( Z\ - T| ) depending on whether we use this upper or lower 
bound for the integral I^. 

By substituting Eq. IV-89, IV-90, IV-91, IV-93 in Eq. IV-88 and by 
expressing the integral 1^ in terms of the approximation given by 

t 

f * 

Eq. IV-97 we obtain for the acceleration pressure drop on the "light" 
fluid region the following expression: 
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where the steady state acceleration pressure drop ([ P* is given by 
Eq. IV-89 . 

IV. 7. 3 The Gravitational Term 

The gravitational term in the momentum equation is given by 
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Substituting Eq. IV-78 and retaining only the first order terms of 
we obtain after integration: 
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The physical meaning of the various terms is : as follows: 

% 

We obtain the steady state gravitational pressure drop by letting 
<£> = in Eq. IV-100, thus in view of the definitions given by 

Eq. IV-36 and IV-75 we have: 
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The second term in Eq„ IV-102 can be expressed by means of Eq. 
IV-102 and Eq. III-22, thus 
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It expresses, therefore, the effect of space lag variation on the 
gravitational pressure drop. 

The third term in Eq. IV-100 can be expressed "by means of Eq. IV-39, 
IV-75 and IV- 103, thus 
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where the mean velocity is defined by Eq. IV-39. This term then 

represents the effect of velocity perturbation in the "light" fluid 
region on the density and therefore on the gravitational pressure drop 
in this region. 

The physical meaning of the last term in Eq. IV-100 can be ex- 
plained again by expressing the integral 1^ in Eq. IV-101 by its upper 
and lower bounds (see Appendix C) 
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which in view of Eq. IV-39, IV-75 and IV-103, can be expressed as 


- 68 - 


si* 

\>u 


rv -5fo-X,) 

e cfU, < T < JL 


A 


-sz b 


s- 1 ^ U, 


ti € su, 


M, 


IV-106 


A simple expression for the integral can be obtained by using the same 
approximation which was used in deriving Eq. IV~97. Thus, if we let 

M ( I in Eq, IV-lOlwe obtain after integration the following 


approximation 
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comparing Eq. IV=107 with Eq. IV“106 it can be seen that has a value 
which falls between the two bounds given by Eq. IV-106. The last term in 
Eq. IV-100 expresses therefore, the effect of the inlet velocity perturbation 
on one density (see Eq, IV-79) and on the gravitational pressure drop in 
the "light” fluid region. Furthermore, this effect is delayed by a time 
delay equal to or depending on whether we use the upper or lower bound 
for the integral I^. 

By substituting Eq. IV-102, IV-103, IV-104 in Eq. IV-100 and by 

expressing the integral I 4 in terms of the intermediate approximation 
given by Eq. IV-107 we obtain for the gravitational pressure drop 
in the "light” fluid region the following expression: 
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where the steady state gravitational presure drop is given by Eq. 
IV-102. 

IV. 7. 4 The Frictional Pressure Drop 

The frictional term in the momentum equation is given by 
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Substituting Eq. IV-78 and I.V-27 in Eq. IV-109 and retaining only the 
first order terms in ^ we obtain after integration the following ex- 


pression 
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where the integral 1^ is given by 
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The physical meaning of the various terms is as follows . 

We obtain the steady state frictional pressure drop by letting 
^ = 0 in Eq. IV-110, thus in view of Eq. IV-35 we can write 
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The second term in Eq. IV-110 can be expressed by means of Eq. IV-112 
and Eq. Ill -22 thus: 
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It expresses therefore the effect of the space lag variation on the 
frictional pressure drop. 

The third term in Eq. Tv-110 can be expressed in terms of Eq. IV-^O 
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This term expresses therefore the effect of velocity perturbation in the 
"light" fluid region on the frictional pressure drop. 

Similarly the fourth term in Eq. IV-110 can be expressed as 
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In view of Eq. IV-79 this term shows the effect of the velocity 
perturbation in the "light" fluid region on the density perturbation and 
therefore on the frictional pressure drop in this region. 

The physical meaning of the last term in Eq. IV-110 can be explained 
again by expressing the integral 1^ in Eq. IV-111 by its upper and lower 
bound (see Appendix C) thus 
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which in view of Eq. IV-115 and IV-35 can be expressed as: 


ji 

5-^ <U,> 



-siT*-r.j 

e 


St 

s-n 


A? " Srt n 

^ € cTu 

<V>,> 


IV-117 


A simple expression for the integral can be also obtained by using the 
same approximation that was used in deriving Eq. IV-97 and Eq. IV-107. 
Thus, if we let in Eq. IV-111 we obtain after integration 
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The last term in Eq. IV-110 expresses therefore the effect of the inlet 
velocity perturbation on the density (see Eq. IV-79) and therefore on the 
frictional pressure drop in the ''light" fluid region. Furthermore this 
effect is delayed by a delay time equal to or (T^*-£ ) depending on 

whether we use the upper or lower bound for the integral. 

By substituting Eq. IV-112, IV-113, IV- 114, IV-115 in Eq. IV-110 and 
expressing the integral 1^ in terms of the approximation given by 
Eq. IV-118 we obtain for the frictional pressure drop in the "light" 
fluid region the following expression: 
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where the steady state frictional pressure drop ^ 11 is given by 
Eq. IV-112. 

IV . 7 . 5 The Exit Pressure Drop 

We can include the effect of the exit pressure drop in the momentum 
equation. For this purpose we shall define by k., ? the coefficient for the 
exit losses , then the exit pressure drop can be expressed as 
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By substituting Eq. IV-69 and Eq. I.V-27, both evaluated at ^ and 

by retaining only the first power in we obtain: 
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We obtain the steady state exit pressure drop by letting £ 
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in Eq. IV-121 thus 
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Consequently Eq. IV-121 can be expressed as 
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The second term in Eq. IV-123 represents the effect of velocity 
perturbation in the "light" fluid region on the exit pressure drop. 
The last two terms in Eq. IV-122 can be expressed as 
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where we have taken into account Eq. IV-81. Consequently the last two 


terms express the effect of the density perturbation on the exit pressure 
drop . 

IV . 7 . 6 The Integrated Momentum Equation 

By adding Eq. IV-86_, IV-98, IV-108, IV-119 and IV-123 we obtain the 
integrated momentum equation for the light fluid region thus 
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By adding the momentum equations for the "heavy" and "light" fluids 
we shall obtain the momentum equation for the system whence the characteristic 
equation for predicting the onset of unstable flow. This will be done 
in the chapter that follows. 

IV. 8 Comparison With Previous Results 

Before we proceed with the derivations of the characteristic 
equation, it is of interest to compare the results derived in this chapter 
with those reported previously in (49 , 50, 51, 53, 55}. In this section 

we shall make comparison with the results of (49, 50 and 51] whereas in 

* 

the section that follows we shall compare the present results to those 
of (53, 55}. 

It was already discussed in Section 1.3 that the assumptions made 
in the present analysis as well as the general formulation of the problem 
are the same as those reported previously the Wallis and Heasley (50} and 
Bourd (51} for boiling, two phase system. It was also noted in Section 
1.3 that the present analysis differs from those reported in (49, 50 and 
51} in the following respect: 1) the constitutive equation of state is 
different and 2) the characteristic equation is different. 

The analyses of (49 , 50 and 51} were derived for boiling systems, 
the present investigation is applicable to both subcritical and super- 
critical pressures. It is emphasized here again that neither this in- 
vestigation nor those reported in (49, 50 and 51} take into account the 
effect of relative velocity between the two phases in the boiling region 
at subcritical pressure.* If the effects of the relative velocity are 

*The conditions under which the effects of relative velocity can be 
neglected are discussed in more detail in (551. 


to be taken into account then the momentum and the energy equation, i.e., 
Eq. IV-2 and IV-3 must be modified. Furthermore, a diffusion equation 


should be added to the field equations describing the process. An in- 
vestigation along these lines will be reported separately. 

If, in the b:oiling region, we express the reaction frequency _TL by 
means of Eq. IV-23, then the density given by Eq. IV-65 becomes identical 
to that derived first in (49} and to those in (50, 51, 55} using different 
approaches. We shall examine now Eq. IV-66 which can be expressed also 
as : 

W) _ _ JL. e ~ T V srt & 

c; rt , 5--a ?(. a, (JI s-sl ^ iv 

IV-126 


whence, in view of Eq. IV-65, and IV-30, the perturbation can be written 
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By adding and subtracting £ / & we can express this relation as 
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If we replace now by Eq. IV-23, then Eq. IV-128 becomes identical 


to Eq„ 38 in the paper by Wallis and Healsey (50) derived using a different 
approach. 


We can insert Eq, IV-65 in Eq. 1V-126 and express the latter as 
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if the density terms which appear on the right hand side of Eq. IV-129 
are .approximated by the steady state relations i.e, ? by 
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as was done in (51) we obtain 


gW) u 
t 


(- 


^lii 


-t 


Sin / -S^ ✓ _ 

££ JI ' +Jl ^ / Ml 


1- 


Ut 


5(S-si) 




i. _ -5. 

e / u, \* 




<>-ri 




IV - 1 3 1 


ich is equivalent to Eq. 5* Appendix A of Boure’s report (51). 

Apart from the difference in the equations of state used in this 
analysis, the difference between the present results and those of (50 * 51) 
is in the handling the momentum equation. In (50) the momentum equation 
was not integrated along the duct, it was first integrated by Bour6 (51). 
Indeed , it can be shown , that after some rearrangement, Eq. IV-88, IV -100 and 
IV-II0 can be put in the form of those given in (51). In (51) the integration 
of the momentum equation lead to a characteristic equation in the form of an 
exponential polynomial of the fourth (or higher) order. 
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In the present analysis -we have introduced various definitions for 
the mean, for the average and for the log mean density as well as for the v@l~ 
©cities in the "light" fluid region which enabled us to give physical 
interpretation to the various terms in the integrated momentum equation. 

It will be seen in what follows that these relations, together with the 
approximation used in deriving Eq. IV-97, IV-107 and IV-118^ result in a 
characteristic equation given by an exponential polynomial of the third 
order. It will be seen also in what follows that these results will 
enable us to derive stability criteria and stability maps which, 
previously, were not available in the literature. 

IV. 9 The Density Propagation Equation 

It is of interest to note an alternate way for determining the density 
perturbation. 

If we substitute Eq. IV-21 in Eq. IV-19 we obtain: 

This equation was called the energy equation in (51) where it was first 
derived. Several remarks are relevant here. 

We note first that Eq. IV-132 predicts the propagation of the 
density caused by the source term it- A " void propagation equation" 
was formulated in (53 and 55) in terms of kinematic waves which predicts 
the propagation of density perturbations through a two-phase system. 
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This void propagation equation takes into account the effect of the 
relative velocity between the two phase as well as the effect of the 
non-uniform velocity and concentration profiles in the two phase mixture. 

It can be easily shown that if these effects are neglected the void 
propagation equation can be reduced to Eq. IV-132. 

We note also that Eq. IV-132 is of the same form as the continuity 
for a given species in a multicomponent , chemical reaction system. In 
chemical kinetics the source term in Eq. IV-132 is referred to as the 
reaction frequency. It is for this reason that in £53, 55) the term 
was called the "characteristic frequency." 

Finally, we note that Eq. IV-132 is a first order partial differential 
equation which can be solved by the standard method used in Sections III-3 
and IV-4. Indeed following this procedure, used already in f53 and 55), 
one can derive Eq. IV-66 and Eq. IV-68. 
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V. The Characteristic Equation 


V. 1 . The Momentum Equation for the System 

The momentum for the "heavy" fluid is given by Eq. III-33, whereas 
that, for the "light" fluid is given by Eq. IV-125. By adding these two 
equations j we obtain the momentum equation for the system. 

We note that if the downstream pressure is constant we can ex- 
press the overall pressure drop, i.e., the external pressure drop of the 
system as a steady state term and a pressure perturbation caused by the 
inlet flow. Thus 
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where the aecond term on the right hand side is determined by the pump 
characteristics and has a negative value. 

By adding Eq. III-33, Eq. IV-125 and Eq. V-l we obtain the integrated 
momentum equation for the heated duct, thus 


~ A Po^ 4- -f- A 4— 4— t 4- A — f- ^ ^3 V *T 

V-2 



t 


-a 


S-SL 


Af« 


u 


tv* 


A Pi 




H 






-h 


^ 


ju« _ 


/l 


vS-XL 




'M, 


ai, 




&* i > 
Mi 


$y 


-X^-ZO 

« cfu, 


Uli 


i'-* 


<T A ■=. o 


We obtain the steady state pressure drop for the system by letting 
the perturbations go "to -zero, thus 

ilfejc -= Afoj aPu -f- A A 

V-3 


By subtracting Eq. V-3 from Eq. V-2 we obtain the perturbed form of 
the momentum equation, thus 
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The formulation is now essentially complete because Eq. V-4 is the 
expression which gives the response of the system to the initial flow 
pertubation as function of the influence coefficients defined below. 

The influence coefficients F 1 and F 2 represent the mass of the "heavy" 
fluid and of the "light" fluid respectively , thus 




f\ ? ■ ily 
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The coefficient F3 describes the effect of the inlet flow variation 
on the pressure drops in the "heavy" fluid region, thus 



This coefficient, which is well known from studies of the transient 
response of single phase flow systems, has always a positive value. 

The coefficient F^ shows the effect of the velocity perturbation in 
the "light" fluid region on the pressure drops in that region, thus 


Aft, 2 Afij 2 .^,, V-8 

+ 4 - — — 

iCUj) w 3 

It is of considerable importance to note that each pressure drop is 
differentiated and is weighed therefore by a different velocity. This 
important result is a consequence of the integration of the momentum 
equation, i.e., of the distributed parameter analysis. We note that in 
the "lumped" parameter analysis the three pressure drops in Eq. V-8 would 
have been divided by ;the same velocity, say by the velocity at the 

exit from the test section as is most often the case for analyses reported 
in the literature. 

The influence coef f icients F^ and F^ are given by 
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It can be seen from Eq. XV-69 and Eq. V-2 that these two coefficients account 
for the effect of the density perturbation on the various pressure drops in 
the "light 1 ' fluid region. Note, that the density perturbation depends on both 
V) and M ( . Two observations are noteworthy. First, the coefficient F^ 
shows that the effects of the velocity perturbation on the "light" fluid 
region' are weighed by various velocities. This, again, is a consequence of 
the distributed parameter approach. Two, the exponential which multiplies 
the coefficient Fg indicates that the effects of the inlet perturbation are 
delayed by the delay time . 

Finally the coefficient F^, defined by 


j- | [ 7) 

1 _a 1 (1-x) 


a 


mH) 


v-u 


shows the effect of the space lag perturbation on the acceleration pressure 
drop in the "light" fluid region. It is important to notice here that In 
Eq. 111=33 and Eq. IV“125 all other terms which are differentiated with 
respect to the length cancel each other in the momentum equation for 


the system. This result could not have been anticipated in a "lumped" 
parameter analysis. Indeed^ in several studies of boiling systems using the 
"lumped" parameter approach these terms were introduced and retained in 
the analysis. In view of the foregoing^ the results and conclusions based 
on such formulations can be considered as spurious. 

By introducing Eq. V-5 through V-ll in Eq. V-4 the perturbed momentum 
equation for the heated duct can be expressed by 
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Before deriving the characteristic equation it will be instructive 
to express the perturbations in Eq. V-12 in terms of the perturbations of 
the inlet flow and of the space lag. Taking into account Eq. IV-30 we can 
express Eq. V-12 as 




d ft > 


d t 


+ 


^ F}> F^j S^i 


si. 


a - su 




\FF-F b e - -\Sk, 


V-13 




d fk 
d h 


-{ 


Fl 


si 


S -JV 


fr-F F n ^siS^ = O 


It can be clearly seen from Eq. V-13 that the dynamic response of the 
heated channel depends upon both the inlet flow perturbation and the 
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variation of the space lag. This latter effect is an example of a 
fluctuation which occurs inside the system. It was discussed in Sections 


II-4 and II-7 that such fluctuation have a destabilizing effect. The 
destabilizing effect which the space lag variation has in combustion 
systems and in boiling systems has been already demonstrated in £48, 52) 
and £50, 51) among others. Equation V-4 shows that, at supercritical 
pressures, the space lag variation has a similarly destabilizing effect. 
•Furthermore, the negative sign in the third term on the left hand side of 
Eq. V-13 shows the destabilizing effect of the inlet velocity perturbation. 
We have noted already that this effect stems from the density perturbation 
in the "light" fluid region. 

V. 2 The Characteristic Equation 

In view of the definitions of the inlet velocity perturbation and' of 
the space lag perturbation given by Eq. III-7 and Eq. IV-30 respectively, 
we can express Eq. V-13 as 
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From this relation we obtain the characteristic equation by noting that 
since £ ^*0 the sum of the terms within the bracket must be equal to 
zero. Thus, after multiplying by (S-JL) and after some rearrangement we 


obtain the characteristic equation for the heated duct: 
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It can be seen that the characteristic equation is a third order polynomial 
with two time delays. From the definitions of the influence coefficients 
we have the following relations for the various terms which appear in 
Eq. V-15. 
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It was discussed in Section II-6 that the characteristic equation 
predicts the value of S as function of the pressure terms given by Eq. 

V-16 through V-22. In general s is a complex number & = & -rLCO f the real 
part gives the amplification coefficient of the particular oscillation 
mode, whereas the imaginary part represents the angular frequency OS . 

Since the original perturbation of the inlet velocity was assumed to be 

ci , str 

of the form J a given oscillating mode will be stable, metastable 

or unstable depending on whether the real part of S is less, equal or larger 
than zero, i.e., whether or a>o. 

A general study of the flow behaviour entails an investigation of 
conditions leading to aperiodic as well as to periodic phenomena. The 
first pertains to the possibility of flow excursion whereas the second 
pertains to the onset of flow oscillations. Following the standard pro- 
cedure we shall study aperiodic phenomena by considering the case of 
S = a with 6j = 0. Again, following the standard procedure we shall study 

i 

periodic phenomena by setting (a = 0^ O) £ o ) in the characteristic 

equation. Such an approach will enable us to determine the stability 
boundary which defines regions of stable and of oscillating behaviour in 
a stability map. In the study of the oscillatory phenomena we shall con- 
sider separately the case of high subcooling and the case of low subcooling. 
The stability problem at intermediate subcoolings will be considered in a 
separate report. 
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VI . Excursive Instability 


VI . 1 Derivation of the Stability Criterion 

The study of excursive, i.e.,, of aperiodic instabilities is con- 
ducted by considering the exponent $ of the velocity perturbation to be 
real, i.e., by letting the angular frequency 60 of the disturbance be zero. 
It follows then from Eq. V-15 that for small values of S , we have the 
following relation: 
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whence after rearrangement: 
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Vl-3 


Equation VI-2 predicts the value of the exponent S in terms of the 
influence coefficients. Since the inlet velocity perturbation is of the 
form of , and since the coefficient A* is positive and exponent 

S is real, Equation VI-3 indicates that the flow will be stable, i.e., 
the disturbance will decrease with time if B* is positive, thus from 
Equation VI-2 

Fj, + ( i— - F r ( 1 -JttL) + fj, -7> S). % > o 

VI -4 

If B is negative then Equation VI-3 indicates that s will be real 
and positive, consequently any flow disturbance will be amplified 
with time resulting in flow excursions. Substituting the definitions 
for the influence coefficients given by Equation V-5 through Equation 
V- 11 we can express Equation VI -4 in terms of steady state pressure 
drops, thus 



This inequality can be cast in a compact form by means of the 
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identities listed below: 
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These relations can be easily derived from the definitions of the steady 
state pressure drops* 

- ; « 9 3 - f i; ; : 


Substituting Eq. VI-6 through Vl-12 in Eq. Vl-5 we obtain the stability 


criterion: 
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which can be expressed also in terms of the total mass flow rate W, thus 
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For boiling systems , this simple criterion was first derived by 
Ledinegg C24J using a different approach, it was analyzed further in 
£25 through 47j and t 51 3- The results of this analysis show that this 
"Ledinegg instability" can occur also at supercritical pressures. The 
significance of the stability criterion given by Eq. Vl-14, can be best 
analyzed by considering the steady state relation fpr the heated 

duct. This will be done in the section that follows. 

VI . 2 Significance of the Stability Criterion 

If, for simplicity, we neglect the effect of the gravitational force 
and if we express the steady state pressure drops in the heated, duct in 

t 

terms of the total mass flow W, and of the total heat input 0 , we have 

the following relations: 
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The total pressure drop for the heated duct is obtained by adding 
Eq. Vl-15 through Vl-19^ thus 
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where the coefficients a, b and c are given by 
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Vl-22 



It should be noted that Eq. Vl-20 is applicable to subcritical as 
well as to supercritical pressures „ By assigning the proper expression 
to (dW* 1 ), which we obtain from the equation of state, we can 
differentiate the process of boiling at subcritical pressures from the 
process of heat transfer at supercritical pressures. Thus, for boiling 
at subcritical pressures we have from Eq. IV-15 
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whereas at supercritical pressures we obtain from Eq. IV-8 
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When Eq. Vl-24 is substituted in Eq. Vl-21, Vl-22 and Vl-23, then Eq. Vl-20 
becomes the pressure drop relation first derived and discussed by 
Schnackenberg (25) and Ledinegg (24) for boiling systems. For super- 
critical pressures Eq. Vl-20 was derived by the writer(63) (see also 
Appendix B) . 


It can be seen from Eq. Vl-20 that whether in boiling at subcritical 
pressures or in heating at supercritical pressures the steady state 
pressure drop in the heated duct has the same cubic dependence upon the 
total mass flow rate. This important conclusion from analysis is indeed 
supported by the experimental data reported by Krasiakova and Glusker (18} 
for water in forced flow through a circular heated duct. Figure Vl-1, 
which is reproduced from Ci8), shows that in boiling at subcritical 
pressures (P = 140 bars) as well as in heating at supercritical pressures 
(P = 226 bars) the pressure drop in the heated duct has the same cubic 
dependence upon the mass flow rate. It could be anticipated therefore that 
the system will have similar dynamic characteristics at these two pressure 
levels. This is indeed the case as it will be shown later. 

The significance of the stability criterion given by Eq. Vl-J-4 can be 
best analyzed by plotting Eq. Vl-20 together with the pump characteristic 
on the same graph. Figure Vl-2 shows such a plot together with three 
possible flow delivery characteristics , i.e., 1) constant pressure drop 
delivery system,, 2) constant flow rate delivery system and 3) delivery 
system specified by the pump characteristics. The intersection of the 
pressure drop for the heated duct with the pressure drop curve of the 

delivery system determines the operating point of the system. The 

' ’■ - •• ■ - «■ ■"'•••' • 

stability criterion given by Eq. Vl-14 indicates that for some of these 

operating points the system, may be unstable with respect to some small 
flow perturbations. In order to show this we shall consider each flow 
delivery system separately. 
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FIGURE 3ZE-I HYDRAULIC CHARACTERISTICS OF WATER FLOWING THROUGH A HEATED 
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FIGURE 3ZT-2 EXCURSIVE FLOW INSTABILITY 



VI. 2.1 Constant Pressure Drop Supply System 


The operating point for a constant pressure drop delivery system are 
indicated by points 1, 2 and 3 in Figure Vl-2. The stability criterion 
given by Eq. Vl-14 indicates that operation at points 1 and 3 will be 
stable whereas that at point 2 will be unstable. For example, if a points 
1 and 3 the flow is slightly increased the pressure drop of the heated 
duct increases, i.e., the "demand" curve of the system increases above the 
"supply" curve of the delivery, consequently the flow will return to its 
original value. Similarly, if at points 1 and 3 the flow is decreased 
the pressure drop of the delivery will be above that required by the 
heated duct resulting in an increased flow and return to the original 
operating point. However, the operation at point 2 will be unstable with 
respect to either a flow increase or a flow decrease. If the flow is 
slightly increased at point 2 the external system supplies more pressure 
drop than that required to maintain the flow. Consequently the flow rate 
will increase until the new operating point is reached. Similarly, if the 
flow is decreased at point 2 more pressure drop is required to maintain 
the flow than is being supplied by the delivery system. Consequently 
the flow will decrease until the new operating point 3 is reached. 

The preceeding considerations can be expressed in a mathematical form 
by noting that for a constant pressure drop delivery system Eq. Vl-20 
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It can be seen from Eq. Vl-26 that flow stability requires an in- 
creasing pressure drop with flow rate. This is indeed the characteristic 
of most flow systems. However , the negative term in Eq. Vl-27 indicates 
that for boiling systems as well as for systems at supercritical pressures 
the pressure drop may decrease with flow rate resulting in flow excursion. 

Instead of the stability criterion given by Eq. Vl-26 one can introduce 
the coefficient of stability S, apparently first proposed by Schnackenberg £25J 
and defined by 
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which in view of Eq. Vl-20 and Vl-27 can be expressed as 
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where the coefficients a, b and c are given by Eq. Vl-21, Vl-22 and Vl-23. 

As observed by Schneckenberg £25J the stability coefficient defined 
by Eq. Vl-28^ represents the per cent change in the pressure drop by a 1% 



- 99 - 


variation in the mass flow rate. It can be seen from Eq. VI-28 and VI-26 that 
for stable flow S must be positive 3 thus 


S> o 


VI-30 


Vl.2.2 Constant Flow Delivery System 

The operating point for a constant flow delivery system is given by the 
intersection of the pressure supply with pressure demand curves. It can be 
seen from Figure VI-2 that for such a system 
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whence Eq. VI-14 indicates that for such a system no flow excursions are possible. 
VI. 2. 3 Delivery Specified by Pump Characteristics 

The operating points for a system whose flow delivery is specified by 
the characteristics of the pump are shown as points 4, 5 and 6 on Figure VI-2. 
Using exactly the same arguments as those used in discussing a constant pressure 
drop delivery system a it can be shown that the operating points 4 and 6 are 
stable whereas operating point 5 is unstable with respect to small flow 
disturbances. At this latter point any flow perturbation will cause a flow 
excursion to either point 4 or to point 6. 


- 100 - 



V3L 8 3: The Effects of Various Parameters and fche Methods for Improving Flow Stability 

The effects which various parameters have on the propensity for flow 
excursions can be evaluated by examing Eq. VI-14, Vl-21, Vl-22, Vl-23 and 
Eq. Vl-27. It can be seen that the variation of any parameter which tends 
to increase the value of the coefficient b given by Eq. Vl-22 will have a 
destabilizing effect. Consequently, increasing the value of the exit pressure 
drop coefficient is destabilizing whereas the flow can be stabilized 

by a high inlet pressure drop, i.e., by appropriate orificing. In view 
of Eq. Vl-24 and Vl-25 it can be also seen that increasing the system 
pressure will have stabilizing effect whereas a decrease in system pressure 
has the opposite effect. Furthermore, the flow can be also stabilized by 
changing the pump characteristics. 

Before closing the discussion of excursive instabilities it will be 
instructive to illustrate the destabilizing effect of the compressibility 
of the fluid in the heated duct. It was discussed in Section 1.3 that the 
instability mechanism which is analyzed in this paper is based on the 
effects of time lag and of density variations in the heated duct. 

For simplicity we shall consider only the effect of the frictional 
pressure drop in a system with zero inlet subcooling, i.e., with A L s* 0 
For such a system Eq. III-20 shows that the space lag is also zero. The 
frictional pressure drop is given by 



where the mean specific volume in the heated duct is obtained from 
Eq. IV-77, thus 
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If we insert in Eq. Vl-32 the expression for the average velocity given 

by Eq. IV-35 and since the space lag is zero, we can express the mean 
specific volume as 
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Since Eq. Vl-32 and Eq. Vl-35 show that both A P and are 
functions of W we can express the stability coefficient defined by 
Eq. Vl-28 as - 
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whence from Eq. Vl-36 we have 
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where ( di ) is given by Eq. IV- 15 or Eq. IV- 8 depending on whether 
we are interested in the subcritical or in the supercritical region* 

It can be seen from Eq. Vl-37 and Eq. Vl-30 that for a system where 
the mean specific volume does not depend on the mass flow rate the flow 
will be stable. For such incompressible flow system the coefficient of 
stability S has a value equal to 2. This is also the maximum value of S 
because when the fricition factor f in Eq. Vl-32 is a function of the 
Reynolds number then Eq. Vl-28 shows that S will have a value less than 
two. For example^ for laminar flow it will have a value equal to unity. 

For a boiling system at subcritical pressures or for a process of 
heating at supercritical pressures Eq. Vl-38 shows that the value of S 
can become negative because of the compressibility of the fluid. For 
such systems Eq., Vl-30 shows that the flow may become unstable. 

In closing it should be emphasized that the density effect per se ; 
can lead to excursive flow instabilities. Oscillatory flow instabilities 
results from a combined effect of time lag and of density variation. This 
will be analyzed in the two chapters that follow. 
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VII . Oscillatory Instability at Low $ubcooling 


VII . 1 The Characteristic Equation and the Stability Map 

In this chapter, and in the following one we shall investigate periodic, 
i.e,, oscillatory flow phenomena. For this purpose we shall assume that the 
exponent <S of the inlet velocity perturbation is given by & = L0J where the 
angular frequency CJ , is a root of the characteristic equation, i.e., of Eq. V-15. 
In this chapter we shall consider the case of low subcooling, whereas, in the one 
that follows we shall consider the case of high subcooling. 

For the case of low subcooling the characteristic equation, i.e., Eq. V-15, 
can be simplified by recalling that for low subcooling the time lag "^b, given by 
Eq. III-19 will be short. Note, that the total transit time , which also 

appears in Eq. V-15 need not be short. This can be seen by considering Eq. IV-63, 
i.e.. 
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Consequently, for short a space lag , and a short time lag Tb , the transit 
time may be long for sufficiently long ducts and/or for low inlet velocities. 
It can be seen from Eq. VII-3 that the effect of time lag will be small if 
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which for subcritical pressures implies 
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whereas, at supercritical pressure this inequality implies 
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When the time lag )\ is short, then in Eq. V-15 the exponential term 
which contains T t , can be expanded and the characteristic equatior reduces 


to 
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This equation can be cast in q dimensionless form by defining a dimensionless 
exponent 
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using this new variable, Eq. VII-7 can be expressed as 
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where the dimensionless coefficients a, b and c are given by: 
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and where the total transit time AT is given by Eq„ VII-7. The coefficients 
a, b and c can be expressed also in terms of the pressure drops, thus 
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Equation VII-9 is a second order exponential polynomial with one time 
delay. Stability maps for such polynomials have been recently presented by 
Bhatt and Hsu Q6.T, 64] . One such map is shown in Figure VII -1, it is in 
the c-b plane with the coefficient "a" as a parameter. The lines for which 
the coefficient "a" is constant are stability boundary curves. For example, 
for given values of the coefficients M b'' and "a", the stable region of variation 
for the coefficient "c" is shown by the line segment AB. The segment CD is another 
stable range for constant values of "b" and of "a". 

Figure VII-1 is the stability map which can be used to differentiate 
the regions of stable operation from the region of unstable, i.e., of oscillatory 
flow in the heated duct. However, because of the complicated nature of the 
coefficients "a", "b M , and "c" which appear on this map, it is rather difficult 
to discuss and analyze the effects of the various parameters. It is desirable, 
therefore, to simplify the characteristic equation in order to obtain simple 
stability criteria. This will be done in the section that follows by neglecting 
the inertia terms in Eq. VII-7. 

VII. 2 Stability Criterion for the Case of Small Inertia 

VII. 2 The Characteristic Equation 

If we neglect the inertia terms F^ and F^ in Eq. VII-7, the 
characteristic equation reduces to its simplest form given by 

-SA Z 
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where the coefficients A and B are given by 
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It is important to note that a characteristic equation of the form of 
a first order exponential polynomial with one time delay describes the onset 
of "chugging" combustion instabilities as shown by Crocco and Cheng £48]j . 

Since Eq. VII-16 is of such a form, we can use the results of Crocco and Cheng [4§J 

.s 

to analyze the flow stability in this problem. The difference between the present 
problem and that of combustion is the physical meaning of the coefficients A and B. 
In this problem they depend on various pressure drops in the system which were 
obtained from the momentum equation. In the combustion problem the coefficients 
are obtained from the continuity equation and depend, among others, on the process 
of combustion. 

We note also that the results of S terming C.62J can be expressed in terms 

of a characteristic equation of the form of a first order exponential polynomial 

with one time delay. However, since Stenning £62J did not formulate his analysis 

rk 

of boiling instabilities in terms of the momentum equation, the coefficients in 
his characteristic equations do not depend upon the pressure drops. 


VII .2. 2 Unconditional Flow Stability 

It was shown by Crocco and Cheng [48} that no matter what the value of 
the time delay AX may be the flow will be unconditionally stable if the co- 
efficients A and B in Eq. VII»16 satisfy the following inequality ^ 
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* The problem was formulated in terms of the continuity and of the energy equation. 
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Because of its importance, we shall define this ratio as the Stability Number N 

s 

In view of Eq„ VII-17 and VII-18, it can be expressed as 


N 


s 


1~ ^ j- ~ F r - flti. ( l~ti ~t Fy - £~sr j 



This stability criterion can be put also in the form of 
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whence upon inserting the values for the influence coefficients in Eq„ VII-27 
we obtain the inequality which must be satisfied for unconditional flow stability 
thus 



VII -22 


This criterion clearly indicates the destabilizing effect of the pressure drops 
in the "light" fluid region and the stabilizing effect of the pressure drops 
in the "heavy" fluid region. 

For once- through systems, when the acceleration and the gravitational 

terms can be neglected, Eq. .VII-22 reduces to 
— 
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If we now approximate ky ^3 we can express Eq. VII-23 as 
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Defining by A ?{. the sum of the pressure drops in the "heavy" fluid region. 
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and by A?^ the sum of the frictional and of the exit pressure drops in the 
"light" fluid region 
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we can express Eq. VII -24 as 
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whence 




A?. 




XI (/-ft) 


Jl(/-ft) 


> 


VII-28 


Inserting now the expressions for the characteristic reaction frequency Si. ? 
for the time lag "Cfe , given by Eq. IV-23 and III-19, respectively, we obtain 
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This inequality can be expressed also in terms of the total mass flow rate and 


of the total heat flow $ . Thus 

2A? f It M J_ t (dtr\[ h Q [j_ 

a (i- —i 1 ) i i - *' ^ 


|\ 

(3 / 




Q (|- 

Q 


7- >1 


VII -30 


Again, we differentiate the process of boiling at subcritical pressure from the 
process of heating at supercritical pressures by using the appropriate equation 
of state, thus at subcritical pressure we use Eq. IV-15, i.e.. 
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whereas at supercritical pressures we use Eq. IV-8, i.e., 
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The implication of Eq. VII-30 will be discussed in Section VII. 3. 
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VII .2. 3 Conditional Stability 

Following again Crocco and Cheng [48~| we can determine the relation 
between the critical transit time and the critical frequencies correspond!! 

to neutral oscillations. Such a relation is obtained by separating the real and 
imaginary parts of Eq. VII-16, thus 
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and 
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where the coefficient given by Eq. VII- 18^ can be expressed in terms of the 
pressure drops thus 
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The stability number N , given by Eq. VII- 20, becomes when expressed in 

s 

*• V • f 

terms of the pressure drops: 
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The critical frequency Wc., is obtained from Equation VII-34 and Equation 


VII -35, thus 
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whereas, the critical transit time A "L is given by Eq. VII-35 and Eq. VII-40, 
thus 



As discussed by Crocco and Cheng, [48^ if the inequality given by 
Eq. VII-19 is not satisfied, then stability is still possible if the angular 
frequency of the perturbation and the transit time satisfy the following 
inequalities 

VII -42 


and 
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The system is intrinsically unstable if the directions of the inequalities 
are reversed. Furthermore, when . 

x T VII -44 

AZ- 


then Eq. VII-16 has an oscillatory solution with an angular frequency C*J c . 


The preceding results can be plotted against the stability number N^, 

b 

given by Eq. VII-20, i.e., by Eq. VII-34. For this purpose we shall define 
also the period of the oscillation by 


T = 


2ir 


CJ, 


VII -45 


We can form now the ratio of the critical transit time to the period 

and express it as function of the stability number N , thus from Eq. VII-45 

s 

and Eq. VII -41 we obtain 




(UO c 

2 . r 


21 ? 


U>L 1 /Y< 


VII -46 


The critical angular frequency can be also expressed as functions of 

N , thus from Eq. VII -40, 
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Similarly, by means of Eq. VII-41 we can express the critical transit 

time as function of N , thus 
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Eq. VII-48, VII-47, and VII-46 are plotted versus the stability number N , in Figure 

. . s 

VII-2. The significance of this map is discussed in the following section. 



VII. 3 Effects of Various Parameters and Methods for Improving Flow Stability 


The effects which various parameters have on the propensity to induce flow 
oscillation at low subcooling can be evaluated by examining Eq, VII-22 or 
Eq 0 VII-30. It can be seen that the variation of any parameter that tends to 
decrease the positive value of the left hand side of these equations has a 
destabilizing effect. For example, increasing the various pressure drop terms 
in the. "light" fluid region has a destabilizing effect. Similarly, an increase 
of subcooling tends to destabilize the flow. Vice versa, an increase of the 
inlet pressure drop or a change of the pump characteristics will stabilize the 
flow. 

Although the preceding results have not yet been tested against experimental 
data, the form of the simplified stability criterion given by Eq. VTI-29, seems 
to be correct. This statement is based on a comparison of Eq. VII-29 with the 
empirical criterion for predicting boiling instabilities recently proposed by 
Serov and Smirnov (66). In the nomenclature of this paper, their criterion is 
given by 
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where a and b are two constants to be determined from experiments, D is the 
diameter of the pipe; V q is the volume occupied by the steam and ( Wi*} jdlP) is 
the variation of the specific volume of the steam with pressure. Consequently, 
the second term on the right hand side of Eq. VII-49 represents the effect of 
compressibility. This effect was neglected in the present analysis. 

It was reported by Serov and Smirnov £66} that Eq. VII-49 was successful 
in correlating data and predicting the onset of flow instabilities in boiling of 
water at pressure of 30, 50, 70 and 100 atmospheres. 


If we neglect the effects of compressibility in Eq. VII-49 and compare it 
to Eq. VII-29 and VII-31, it can be seen that Eq. VII-49 is incorporated in 
Eq. VII-29. We note also that this latter equation is a simplified form 
of Eq. VII-23; i.e. of Eq. VII-9 which are therefore more general and 
complete . 

Further experimental evidence that gives support to the form of 
Eq. IV-29 is shown on Figure VII-3 which is reproduced from the paper 
by Platt and Wood /I 0/ . It can be seen from this figure that either 
increasing the power input and/or decreasing the mass flow rate has a 
destabilizing effect. The same results are predicted by Eq. IV-29. 

Perhaps the result of greatest significance revealed in the present 
investigation is the similarity between the characteristic equations for 
predicting "chugging" combustion oscillations and the characteristic 
equation for predicting low frequency flow oscillations in heated ducts 
at near critical and at super-critical pressures. Since it is well 
known (see for example[48~) ) that "chugging" combustion instabilities can 
be stabilized by an appropriate servo-control mechanism, the results of 
this investigation indicate that low frequency flow oscillation at near 
critical and at supercritical pressures may be also stabilized. This 
important conclusion is demonstrated on Figure VII-2 which shows also the 
effect of various parameters on the propensity toward oscillatory flow. 

It can be seen on Figure VII-2 that even when the stability number 

N s is less than unity, the flow may be stable if the frequency of the 

inlet perturbation is higher than the critical frequency ^ <_ . Similarly, 

the flow can be stable if the total transit time is shorter than the critical 

one. The values of U ^ and of = AT C are plotted in Figure VIII-2 

c 

in terms of the stability number N s and of the coefficient B given by 
Eq. VII-39 Eq. VII -36 respectively. 
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FIGURE HI-3 PERCENTAGE OF HEAT EXCHANGER TEST DATA SAMPLES SHOWING STEADY 


OPERATION VS. HEAT FLUX PER UNIT AREA PER UNIT MASS FLOW RATE, 
REPORTED BY PLATT AND WOOD. 
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The effects which the variations of the various parameters have on 
the flow stability can be evaluated from Figure VII-2 by considering 

I- - Wl ' 

whether the variation results in an increase of the stable region. For 
example, it can be seen from Figure VTI-2 that for a constant value of N s 
an increase of the delay time has a destabilizing effect because for 
sufficiently long delays AL will become larger than We note that 

this quantitative conclusion is in agreement with the qualitative des- 
cription of the destabilizing effect of the time delay presented in Section 
II -4. It can be also seen from Figure VII-2 that increasing the frequency 
of the inlet perturbation at a constant value of N s , has a stabilizing 
effect because for sufficiently high frequency will become larger 
than . Furthermore, Figure VII-2 shows that an unstable flow; i.e., 

a flow for which W4<^ c and AT > can be stabilized by increasing the 
value of the stability number N s . 

We close this section by observing that the foregoing conclusions 
and results are new and have not yet been verified against experimental 
data. If confirmed, then the results of this study provides a method 
whereby stable operation can be insured on an intrinsically unstable 
region. 
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VIII. Oscillatory Instability at High Subcooling 


VIII .1 The Characteristic Equation and the Stability Criterion 

We shall consider now the case of high inlet subcooling which implies a 
long time lag ^4 and a long space lag / . For such system Eq. VII-3 indicates 

that the transit time and the time lag will be of the same order of magnitude. 
Since both time delays are long, we shall neglect the exponential terms in the 
characteristic equation given by Eq. V-15, which reduces then to 
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which can be rearranged and expressed as 
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where the sums _of the influence coefficients are related to the pressure drops 
by the following relations 
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It can be seen that the characteristic equation is a cubic equation of the 
form of 
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where the coefficient a] b, c and d are given by the corresponding terms of 
Eq. VII -2. 

The problem of determining the conditions for neutral stability is solved 
again by substituting S = t UJ in Eq. VII-7. 


Thus 
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5 



whence upon separating the real and the imaginary parts we have 
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and 
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Consequently for oscillations to be possible the coefficients a,b,c and 
d in Eq. VII-9 and Eq. VII-10 must satisfy the following relation: 
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whence, the values of the influence coefficients must be such as to satisfy the 
following expression: 
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It can be seen from Eq. VIII-4 and Eq . VIII-6 that, unless the effects 
of inertia or of gravity become dominant, the right hand side of Eq. VIII-12 
is a positive quantity. Consequently, Eq. VIII-12 indicates that 
oscillation can occur only if 
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Therefore} the flow will be stable if 
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In view of Eq. VIII-5, this inequality can be expressed also as; 
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For oscillatory flow, Eq. VIII-13 and Eq. VIII-9 indicate that the angular 
frequency will be given by 
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which, when expressed in terms of the influence coefficients, becomes 
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It should be noted, again, that the values of these influence coefficients should 
satisfy Eq. VIII-16, i.e. Eq. VIII-12. 
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-VIII. 2 Effects of Various Parameters and Methods for Improving Flow Stability 

The effects of the various parameters can be evaluated by examining the 
inequality given by Eq. VIII-15. It can be seen that any variation which tends 
to increase the value of the left hand side of this equation will have a stabilizing 
effect. Thus, the flow can be stabilized by increasing the pressure drops in 
the "heavy" fluid region, whereas it will be destabilized by increasing the 
pressure drops in the "light? 1 fluid region. 

The effect of subcooling can be evaluated by comparing Eq . VIII-14 and 
Eq. VXEI-15 with Eq. VII-20 and Eq. VII-39. Since the velocities in the "light" 
fluid region are higher than the inlet velocity it can be seen from such a 

comparison that the inequality applicable at high subcoolings, i.e. Eq. VIII-14 

I 

is less restrictive than that corresponding to low subcoolings, i.e., than 
Eq. VII-20. Consequently, the flow is more stable at high subcoolings. 

However, since Eq. VII-20 indicates also that an increase in subcooling destabilizes 
the flow, we conclude that this destabilizing effect must go through a maximum 
at intermediate subcoolings. For boiling systems, this conclusion is in agree- 
ment with the experimental results of Gouse (67) who was apparently the first 
to notice this effect. At super critical pressures, experimental data, which 
could be used to test this conclusion, are not yet available. 


IX. DISCUSSION 


The instability mechanism investigated in this paper was based on the 
destabilizing effects of time lags and of density variations in the heated 
duct.* It was shown that, in the near critical and in the supercritical 
region, these destabilizing effects can induce flow excursions as well as 
flow oscillations. 

The characteristic equation, i.e., Eq. V-15, which predicts the onset 
of these instabilities is given by a third order exponential polynomial 
with two time delays. Because of its complex nature this equation was not 
solved at this time. Instead, simplified stability criteria were sought 
and derived by assuming that the inlet subcooling was either low or high. 
This approach seemed preferable for several reasons. 

First, the simple stability criteria are more instructive and helpful 
for gaining an understanding of the essential nature of the instability. 

Two, the result shows that the dominance of a particular parameter re- 
sults in a particular angular frequency of oscillations (see Eq. VII-40 and 
VIII-17). Consequently, the cause of instability can be determined from a 
trace of the flow oscillation. 


*Other mechanisms which may induce flow oscillation were discussed in 
Section II-7. 
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Finally, simplified stability criteria such as Eq. VI-20, VII-22, 

VIX-42 and VIII-15 are more amenable to a qualitative study of the effects 
which variations of the various parameters may have on inducing or on pre- 
venting flow excursions and/or flow oscillations. Indeed, only if the 
results from Such a study are in agreement with experimental observations, 
a detailed quantitative solution of the more complicated characteristic 
equation can be justified. 

It was discussed in Sections VI-2, VII-3 and VIII-2 that the pre- 
dictions based on the simplified stability criteria are indeed in qualitative 
agreement with the experimental data. This agreement warrants therefore a 
more complete study of the characteristic equation together with a quantita- 
tive comparison with the experimental data. 

Last but not least the simple criteria are most useful in indicating 
the improvements and changes in the design or in the operation of the system 
which would insure stable flow. Several such improvements were discussed 
in Sections VI-2, VII-3 and VIII-2. It was noted there that the results of 
this study indicate that low frequency thermally induced flow oscillations 
in the near critical and in the supercritical pressure region, could be 
stabilized by an appropriate servo-control mechanism. Whether this important 
conclusion is indeed correct remains to be shown by future experiments. 
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Appendix A 


The Near— Critical Thermodynamic Region 
The success of an investigation concerned with predicting or in- 
terpreting the behaviour of a thermo-dydraulic system depends on the 
availability and on the accuracy of data giving the values of thermodynamic 

and transport properties of the fluid in the region of interest. It is 

/ 

the purpose of this appendix to summarize, briefly, the status of present 

understanding of thermodynamic phenomena that take place in a region near 

the critical thermodynamic point. For additional discussion, the reader 

is referred to the extensive reviews by Rice (Al) and by Hammell (A2) . 

Consider a fluid at a pressure slightly above the critical pressure 

flowing through a heat exchanger. If the temperature of the fluid at the 

entrance is considerably below the critical temperature, i.e., T « T , the 

c 

fluid will have a density close to that of a liquid whereas at the exit, 

if the fluid temperature is considerably above T , the density will ap- 

c 

proximate that of a perfect gas. Consequently, in passing through the 
heat exchanger the fluid will undergo a change of properties from a. liquid- 
like fluid at the entrance to a gas-like fluid at the exit. Since the 
properties of the fluid will affect the performance of the system is 
becomes necessary first n to examine the nature of this change and then to 
express it quantitatively. 

At subcritical pressures the presence of two phases is distinguished 
by a difference in density and by the existence of an interface between 
the phases. At supercritical pressures such a distinction cannot be made 
because at these pressures as well as at the critical one the interface, 
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the heat of vaporization., as well as the surface energy, all vanish. 

There is no general agreement as to the structure of the medium and 
of the mechanism of phase transition in the critical and in the supercritical 
region. Different explanations and descriptions are advanced by different 
authors . 

Some authors like Rosen (A3) and Semenchenko (A4) analyze the thermo- 
dynamic characteristics of a medium in the supercritical region by assuming 
an equation of state like the? Van der Waals' or the Dieterici equations. 

Kirschf elder , Curtis and Bird (A5) describe the fluid in the neighbor- 
hood of the critical point as consisting of a large number of clusters of 
molecules of various sizes. The system can be idealized by assuming that 
the density can be described by a distribution function which has for its 
two limits the densities of the two phases. The fltictuation in density, 
which can be expected from the theory of fluctuations, becomes very large 
in the vicinity of the critical point. These large fluctuations and the 
formation of molecular clusters in the neighborhood of this point result 
in a large increase of the specific heat at constant volume. 

Mayer and co-workers (A6) propose a theory of condensation based on 
the cluster theory of imperfect gases from which they predict the existence 
of an anomalous region above the temperature of the usually observed critical 
point. This region extends up to the highest isotherm for which (1) V /l)v) T , 
has anywhere a zero value. In this region, isotherms exist having no vari- 
ation in pressure over a finite density range, but having at all densities 
continuous derivatives with respect to pressure. Various aspects of this 
theory are discussed further in (A5) . 
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A great number of authors distinguish two phases in the supercritical 
region: a heavy, liquid-like phase and a light, gas-like phase. The 

difference between their results stems from the different approaches used 
to locate the boundary between the two phases and from the different 
descriptions of the characteristics of the phase transition. 

In a preceeding section we have discussed already Goldman's (A7-A8) 
descriptions of the supercritical region and of the similarity between the 
heat transfer and flow processes at supercritical pressure and those that 
take place at subcritical pressure during the process of boiling. However . 
Goldman did not formulate, quantitatively, the problem nor did he say how 
and where to locate the boundary or the region between the liquid- like and 
the gas- like phase. 

Following Goldman, Hendricks et al (A9) consider "boiling-like" 
phenomena at supercritical pressures and, in analogy with boiling, they 
introduce a specific volume for the fluid of the form of Eq. Al. 


^ = or + 


x 




(A- 1 ) 


In place of the quality they introduce a weighting function for the heavy 
and light species. However, no reference is made in their paper as to how 
to determine, quantitatively, this distribution function. 

In numerous textbooks (A- 10) among others, the boundary between the 
liquid and the gas in the supercritical region is taken to be the critical 
isotherm. Other authors like Thiesen (A-ll), Trautz and Ader (A-12) 
among others take the critical isochor for this boundary and consider it as 
the extension of the saturation line into the supercritical region. 
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In the subcritical region various thermodynamic properties such 
as the specific heat, the compressibility, the coefficient of thermal ex- 
pansion and others change discontinuously or reach a maximum value at the 


coexistence, i.e., the * aturation line. This line can be therefore looked 
upon as the locus of points for these discontinuities or maxima. Conse- 
quently, numerous authors consider the extension of the saturation line 
into the supercritical region to be the line which is the locus of points 
where the thermodynamic properties listed below reach a maximum: 



(A-2) 


(A- 3) 


(A- 4) 


(A- 5) 


(A-6) 


Several authors (A-13 - A-17) assume that one single line represents 
the locus of points of all these maxima. This is indeed the case for sub- 
critical pressure where the saturation line is the locus for all discontinuities 
However, the experiments of Kaganer (A- 18) and of Sirota and co-workers (A-19) 



show that this is not the case but that for a given supercritical pressure 
different thermodynamic properties reach a maximum value at different 
temperatures. Thus, for each of the thermodynamic properties, i.e. , 
specific heat c^, the coefficient of thermal expansion, etc. , there is a 
different line which represents the locus of the maxima. This raises the 
question which of these lines can be regarded to be the extension of the 
saturation line in the supercritical region, i.e., which of these lines can 
be considered as the boundary between the liquid-like and the gas-like 
phase. 

Plank (A-20) and Semenchenko (A-21) consider the line along which 


25 

D v ‘ 


= 0 


(A- 7) 


to be the extension of the saturation line in the supercritical region. 
Eucken (A- 13), however, takes the curve represented by Eq. A- 2 for this 
extension; whereas numerous authors (A- 8, A- 9, A- 22 - A- 25) take Eq. A- 3. 

Of particular interest to the analysis of this paper are the results 
reported in (A- 14, A-17, A-19 and A-16) which will be therefore discussed 
in more detail. 

Sirota and co-workers (A-19) discuss the transition phenomena at sub- 
critical and supercritical pressures in terms of the Frenkel's theory of 
heterogeneous fluctuations (A-14, A-26). According to this theory in any 
gas at subcritical temperature heterogeneous fluctuations result in the 
formation of molecular complexes which can be regarded as finely dispersed 


- 133 - 


nuclei of a phase within a homogeneous phase. In approaching the saturation 
line the fluctuations increase and "micro-heterogeneities” appear in the 
macroscopic, homogeneous phase. This marks the beginning of the "pre- transi- 
tion region" which is characterized by the fact that various thermodynamic 
properties exhibit variations which become more pronounced as the saturation 
line is approached. This accounts for the anomalous effects of the proper- ' 
ties in the vicinity of the saturation line. At the saturation line the 
properties change in a discontinuous fashion which is a characteristic of 
phase transitions of the first order. As the pressure is increased the 
effect of heterogeneous fluctuations increases whereas the effect of phase 
change, i.e., of the discontinuous change of properties becomes less 
important and disappears at and above the critical point. Since the change 
of phase at subcritical pressure is characterized by an obsorption of energy 
and an expansion of volume the transition at supercritical pressure should 
be characterized by the maximum values of c^ and of the thermal expansion, 
i.e., of ( r Dv/ r "])T)p. See Figures A-l and A-2 which show these properties 
for oxygen at supercritical pressures. However, the authors of (A- 19) show 
from experiments that at a given pressure the two maxima do not occur at 
the same temperature. The values of the maxima for are correlated by 

(A- 8) 
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Figure A1 Specific Heat Versus Temperature for Oxygen at Supercritical Pressures. 
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Figure A2 (3v/0T)p Versus Temperature for Oxygen at Supercritical Pressures. 1 






which is valid for non-polar liquids when p / p cr ^ t 1»5. In the above 

equation, R is the gas constant whereas c is the specific heat for an 

p g 

ideal gas. This equation shows that the value of the maximum c^ decreases 
as the pressure is increased. The temperatures where these maxima occur 
were correlated by 


JL_ = _L_ _ 0.4-87 /o^ JL 

1 1 r ..it ^t.r it 


ic Icrit 

This temperature, denoted here by T 


(A- 9) 


pc' 


is often referred to in the litera- 


ture as either the pseudo-critical temperature or the transposed critical 
temperature o 

• vr “ 

Both Sirota (A-19) and Kaganer (A-18) show that the locus of the maxi- 
mum values of c^ along isobars, i.e,, Eq. A-3, is the extension of the sat- 
uration line in the supercritical region. 

Urbakh (A-17) also considers the effect of heterogeneous fluctuations 
at subcritical and supercritical pressures. He shows that as the temperature 
is increased and the surface tension decreases the heterogeneous fluctuations 
increase and reach a maximum at the critical point. The location of the 
critical point depends on the surface tension; moreover, it can be changed 
by introducing surface active agents. The critical point divides two regions 
which can be distinguished by the nature of the phase transition. At sub- 
critical pressure the transition is characterized by the discontinuities 
of the properties and by the presence of a macroscopic second phase within 
the originally homogeneous phase. At supercritical pressures the second 
phase is finely dispersed in the form of clusters. Furthermore, in this 
region the properties do not change discontinuously but vary in a continuous 
way. At subcritical pressures the effect of heterogeneous fluctuations 
becomes evident in the "pretransition region" as a variation of properties 
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in the vicinity of the saturation line. This is shown in Figure A-3 which 
is the volume- temperature plane for oxygen. At a subcritical pressure, 
say at P^ = 0.9, the line 1' - 2' is the phase transition of the first 
order occurring at a constant temperature. The effect and magnitude of 
the fluctuation in specific volume in the two pre-transition regions is 
shown as the lines 1-1* and 2-2'. The fluctuation 1 - 1* is caused by the 
formation of vapor nuclei in the pre- transition region of the liquid. 
Similarly, 2-2' are the fluctuations caused by the formation of liquid 
nuclei in the pre-transition region of the gas. It can be seen from this 
Figure that at low pressures in the subcritical region the effect of 
fluctuation is negligible when compared to the phase transition of the 
first order. For example, at P = 0.5, they are almost absent. Increas- 
ing the pressure increases the effect of heterogeneous fluctuations which 
reach a maximum at the critical point. At this point and above it^the 
phase transition of the first order vanishes so that only the effect of 
heterogeneous fluctuations remains. Urbakh notes further than with the 
phase transition and the fluctuations are associated energy requirements 
which can be determined from the T - s or v - s diagrams shown on Figures A- 4 
and A-5. At low pressure the only energy required is heat of vaporization 
for the phase transition of the first order, thus 

h f = T ( s 2 ' - s^) (A-10) 

However , as the pressure is increased the energy associated with the 
fluctuation becomes important. At supercritical pressure it is the only 
which remains, and it can be determined either from Figure A-4 or A-5, 
thus 

AL = T ( s 2 - a x ) (A- 11) 
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Frenkel (A- 14* A-26) considers two variations i a transition of the 
first order at subcritical pressure and a transition of the second order at 
supercritical pressure. The first , characterized by discontinuities of 
properties, is described by Glausius-Clapeyron' s equation: 


dP 

dT 


M. 


T (v - v f ) 
o g f 


(A-12) 


and takes place at a constant temperature T . The phase transition of 
the second order takes place over a temperature interval A T = T^ ~T^, 
in which the properties change continuously . In this temperature interval 
both 0 ^ and 0)-v/^T) reach a maximum. FiguresAl andA2 show these varia- 
tions for oxygen at three supercritical pressures. As a generalization 
of the transition of. the first order Frenkel formulates the equivalent 
energy of transition for the second order transition, thus 


AL = T tc Cs 2 ' s l ) = 


A c dT 

^ p 


(A- 13) 


where T^is the temperature corresponding to the peak of c^ and A c^ is 
the value of c^ above the "normal" value, i.e., above the dashed line on 
Figure 1. Similarly, the change of volume is given by 


V 1 = 


£ Vh). 


dT 


(A- 14) 
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Eq. A- 13 and Eq. A- 14 represent the additional increase of volume and the 
additional heat absorbed in going from the liquid-like state to the gas- 
like state at constant pressure. In place of Clausius Clapeyron's equation 
Frenkel uses the equation derived by Ehrenfest (A-27) to describe transitions 
of the second order at the "lambda point" of helium and at the "Curie point" 
of feromagnetic metals, thus 


dP 

dT 


A c 


m 


T tc* 


m 


(A- 15) 


where Ac and A (Iv/^I) are the maximum values of c and of ( v/ T)_ 

^ m ' p p ' P 

above the dashed lines in Figures A1 andA2. Various criticisms which have 
been made with respect to Ehrenfest equation are discussed in (A-28). Also, 
various authors (A-18, A-19) criticize the use of Eq. A-15 for the supercritical 
region because the temperatures where c^ and (D v/ 0 T)p reach their 


respective maximum values are not the same. Consequently, the value of T 


tc 


in Eq. A-15 is somewhat arbitrary. 

Semenchenko (A-4, A-16, A-29) considers the medium in the supercritical 
region to consist of two phases which are separated by a region in which the 
properties change rapidly but continuously. It was already noted that he 
takes the locus of points given by Eq. A-7 to represent the extension of 
the saturation line in the supercritical region. He notes that at subcritical 
pressures the phase transition is accomplished by absorbing an amount of 
energy given by Eq.AlO and by doing an amount of work given by: 


W = P (v 2 - v x ) 




(A-16) 
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However, since in the supercritical region there is no discontinuous change 
of volume and of entropy, Semenchenko notes that Eq. A- 10 and A- 15 must be 
modified and replaced by: 



(A- 17) 


(A- 18) 


For additional discussion of critical phenomena the reader is referred 
to the extensive reviews by Rice (A-l) and by Hammell (A-2) . 

From the preceding review of the present understanding of thermodynamic 
phenomena in the supercritical region we can make the following conclusions: 

1) There is no general agreement as to the structure of the medium 
and of the mechanism of phase transition in the critical and super- 
critical region. 

2) There is a general agreement that large variations of density and 
of specific heat are present. 

3) Most of the authors consider the supercritical region to consist 
of two phases -- a liquid-like and a gas-like phase. 

4) There is no general consensus as., t;o the location of the boundary 
or of the transition region between these two phases, although a 
large number of • investigators consider this demarkation to take 
place along the line which is the locus of points where the 
specific heat at constant pressure reaches a maximum. 
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5) There is no general consensus as to the nature of phase transition 
at supercritical pressures and of the energy required to bring it 
about. Three different methods for evaluating this energy of 
transition have been proposed: 1) the graphical method of 

Urbakh (A-17) resulting in Eq. A-ll; 2) the second order transition 
proposed by Frenkel (A-14, A-26) given by Eq. A-13, and Eq. A-15; and 
3) the pseudo transition region proposed by Semenchenko (A-4, A- 16, 
A-29) given by Ej. A-35 and Eq. A-17. By examining the proposed methods 
and equations, i.e., Eq. A-ll, Eq. A-15 and A-17, it can be seen that 
these different methods will yield different values for the transition 
energy. 

It is evident from the preceding results that the success of any analysis 
concerned with the mechanism of flow oscillations and of heat transfer at 
supercritical pressures will depend to a great extent upon the ability to 
describe more accurately the thermodynamic state of a fluid and the transi- 
tion phenomena that take place at supercritical pressures. 
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Appendix B 


The Steady State Pressure Drop 

In this Appendix we shall derive an expression for the steady state 
pressure drop of a fluid whose properties change fro’ 7 ', a liquid- like at the 
entrance to a gas-like at the exit of the heat exchanger. The derivation 
and the resulting flow excursion criterion applicable to fluids at critical 
and supercritical pressures were first derived by the writer in the Second 
Quarterly Progress Report. They are reproduced here for reasons of completenes 

The pressure drop across a heated length L is the sum of the acceleration, 
pressure drop, the frictional pressure drop and the pressure drops across 
the inlet and exit flow restrictions. Since the pressure drop depends on 
the fluid, it becomes necessary to examine first property changes along the 
heated duct. 

B . 1 The System - Three Region Approximation 

The system analyzed in this Appendix is shown in the Figure B-l. 

6 

A circular duct is uniformly heated at a rate of Q, over a total heated 
length L. Two flow restrictions are located at the entrance and at the exit 
of the heated section. A fluid at an initial temperature T^» i.e., with the 
enthalpy i^, flows at a constant mass flow rate . In passing through the 
heated duct the specific volume and the enthalpy of the fluid increase (See 
Fig.-B-l). The fluid undergoes, therefore, a transformation from a liquid- 
like to a gas-like fluid. 

Figure B-2 shows the V -i relation for oxygen at a reduced pressure of 
P = 1.1. It can be seen from this figure that the increase of specific 



volume from a liquid-like state to a gas-like state occurs gradually over 
an enthalpy interval. 


In order to simplify the problem, we shall assume that the entire 

transformation can be approximated by considerding three regions . In the 

first region of length l f , between stations (l) and (|) in Figure B-l, the 

heavy clusters resemble a liquid. In this region the specific volume of 

the fluid is constant having a value of v^. We shall assume that the com- 

ple te transformation, from heavy to light clusters, takes place within the 

transition length l t , i.e. , between stations (2) and (£) . In this transition 

region the specific volume of the fluid changes from a value of v^ to a value 

of v . The enthalpy change associated with this expansion is given by 
§ 2 

♦ * t 

^ t, , - t T • In the third region of length 1 , the light clusters 
resemble a gas. The specific volume of the fluid in this region can be 
approximated by that of gas and, in particular, by that of a perfect gas. 

It is apparent from the discussion in Appendix A that the initial and 
the final conditions of the transition region, i.e., the conditions at stage ( 2 ) 
and (^) respectively, will depend upon the model selected for describing 
the pseudo-phase transition in the supercritical region. This follows from 
the fact that the temperature or the enthalpy that marks the start of the 
pseudo-phase transition will determine the location of station , whereas 
the location of station will depend on the energy required to complete 
the transition from heavy to the light clusters. In this report we shall 

t 

denote this energy requirement by which can be determined by the best 

three region approximation indicated in Figure B-l. 

As discussed in the preceeding sections, we are considering in this 
report only the effects of density variation on the flow stability. 


Consequently, we shall assume that both the friction factor and the heat trans- 
fer are constant. The first assumption is quite reasonable if the flow remains 
turbulent throughout the duct. The limitation of the second assumption may 
become significant if variations of transport properties in the transition 
region have an important effect on the stability. We note, however, that 
both assumptions can be removed permitting an extension of the analysis to 
consider the effect of variations, other than density, on the initiation of 
flow oscillations. 

B. 2 The Frictional Pressure Drop 

The frictional pressure drop in the system is given by the sum of the 
frictional pressure drops across the segments 1^, l fc and 1 and the pressure 
drops across the inlet and exit flow restrictions, thus 

+ A\y fw) i-d^cur) i 

For a constant friction factor f, the pressure drop across a segment 
of length 1 is given by 

IT) v, A 1 B-2 

where the lengthwise average specific volume is given by 


i 



Consequently, in order to evaluate the frictional pressure drops for 
the three segments it is necessary to evaluate the specific volume for each 
segment. This can be done by relating first the specific volume to enthalpy 
and then to express the enthalpy in terms of the heated length. This latter 

relation can be obtained from energy considerations. 

% 

Denoting by Q, the total rate of energy addition to the system and by 
Q. the constant heat flux density, we have for a duct 


77 - r * B-4 

where ^ , is the heated perimeter. It follows from Eq. B-4 that 

B-5 


B-6 


Furthermore, for a system with constant mass flow rate 
enthalpy is given by 

W ol i d Q . 

where we have neglected the kinetic energy of the fluid. It follows then 
from Eq. B-7 and B-4 that 

yr <u = b -8 
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the change of 
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where the total length is given by 


L - ^ -t* *+■ ^ 


and in view of Eq. B-5 we obtain 
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Substituting Eq. B-4, B-3 in Eq. B-2, we obtain the pressure drop across 
a heated segment where the enthalpy of the fluid changes from i to i +^i, 
thus 
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, For a three region approximation the relation between v(i) and i 
is shown in Figure B-l. We shall consider now each region separately. 

a) The Liquid-Like Region 

In this region the specific volume of the fluid is constant and equal 
to v^ (See Figure B-l). In the segment of length 1^, the enthalpy of the 
fluid increases from ij- to i^. The frictional pressure drop across 1^ be- 
comes then 
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b) The Transition Region 

In the transition region we shall approximate the relation between the 
specific volume v, and the enthalpy i, by a linear equation. The average 
specific volume in this region can be written then as: 
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Denoting by- 
pressure drop in 


» * ' 

the change in enthalpy, the frictional 
the transitional region then becomes 
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c) The Gas-Like Region 

In view of the assumption that in this region the fluid has the prop- 
erties of a perfect gas we have, for a constant pressure process, the fol- 
lowing expression for the specific volume 


~ H" 


u 


X 

y cj* 




B-14 


Inserting this expression in Eq. B-10 we obtain 
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The change of enthalpy £-^2 ’ can ex P resse d also in terms of the total 
heat input thus from an energy balance 
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Inserting this expression in Eq. B-15, we obtain for the frictional 
pressure drop in the gas-like region the following expression 
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Denoting by a numerical coefficient that takes into account the 
geometry of the restriction and of other losses like vena contracta etc., 
we have the inlet pressure drop 
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Similarly, we define by k^ a numerical coefficient that accounts 
for the geometry and the losses at the exit. The exit pressure can be then 


expressed as: 


K ^ w 


B-19 


which, in view of Eq. B-14 and B-16, can be also written as: 
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B.4 The Acceleration Pressure Droi 


The acceleration pressure drop is given by 



I 


then 
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The total acceleration pressure drop required to accelerate a fluid of 
specific volume at the inlet up to the exit where it attains a specific 
volume v^g is given therefore by 
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Inserting Eq. B-14 in Eq. 24 we obtain 
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or in view of Eq. B-16 we have 


A? « = tr) 1 ^ ^ (' vr ‘i.) ~ 


R / <5 


B-26 


B . 5 The Total Pressure Drop 

The total pressure drop is obtained by summing Eq. B-26, B-20, B-18, 
B-17 and B-ll, thus after some rearrangement 
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where the coefficients a, b and c are given by 
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The form of Eq. E-27 is relevant to the present problem because it shows 
that, for some operating conditions, the pressure drop may decrease with 
increasing mass flow. This consequence of the negative term on the right- 
hand side of Eq. B-27„ 

B. 6 The Two Region Approximation 

Following the derivation of the pressure drop given in the preceding 
section it was observed by Dr. R. Fleming, from the Research and Development 
Center at G.E., that instead of considering a three region approximation as 
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shown in Fig. B-l, the problem can be further simplified by considering a 
two-region approximation indicated in the sketch below 



Two- Region Approximation 


In the two-region approximation the transition region shown in Fig; B-l 

and B-2 is neglected, i.e., 1 = 0, i. = iJ, Sr . = v r . It is assumed, 

t 2 2 w g2 f 

therefore, that the change from a liquid- like to a gas- like fluid occurs in- 
stantaneously in a plane perpendicular to the flow when the enthalpy reaches 
a value of i^ indicated on the sketch above. 

We can further amplify the preceding observation. It can be seen 
from Fig. B-2 that the enthalpy which corresponds to the transition point 
can be approximated by the enthalpy at the transposed critical temperature 
T ; , i.e., by the enthalpy that corresponds to the maximum value of the 
specific heat at constant pressure c^. Consequently, with a two-region 
approximation one can consider that the liquid-like state persists until 
the temperature of the bulk fluid reaches a value that is equal to the 

transposed critical (or pseudo-critical) temperature . Above that temperature 

; -153- : 


the fluid behaves as a gas. Therefore, the transposed critical temperature 

T , can be regarded as the boundary between the liquid- like and the gas- like 
t c 

states. It was discussed already in Appendix A, that both Sirota (A-19) 
and Kaganer (A- 18) have shown that this temperature is the extension of the 
saturation line in the supercritical region. Figure A-l in Appendix A 
shows that the transposed critical temperature increases with increasing 
values of reduced pressure. It can be concluded therefore that the value 
of enthalpy corresponding to this temperature and to the transition point 
shown in Fig. B-2 will also increase with increasing reduced pressures. 

For a two-region approximation the form of Eq. B-27 remains unchanged, 
however, the coefficients a, b and c given by Eq. B-28 and B-29 and B-30 
reduce to: 
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which we obtain by setting = = Oj V / = in Eq. B-28 


B-29 and B-30. 


As noted by Dr. Fleming the use of the two-region approximation simplifies 
considerably the form of the coefficients a, b and c. The three region 
approximation retains however a closer similarity with phenomena that take 
place at subcritical pressure. The transition region shown in Fig. B-l can 
be regarded as corresponding to the boiling region at subcritical pressures. 

The liquid and the gas region in Fig. B-l would then correspond to the pre~ 
heating and to the superheating region in a once- through boiling system 
where the liquid at the entrance is subcooled and the steam at the exit is 
superheated. We have noted already in Appendix A that the enthalpy change 

considered as being equivalent to the heat of vaporization 
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The selection of either the two or three region approximation should 


be determined by the desired simplicity and accuracy. The important result 
is however the fact that, because of the negative term on the right-hand side 
of Eq. 27, there exist a possibility of a decrease in pressure drop with 
increasing flow in the supercritical thermodynamic region. It was shown 
in the body of the report that such a pressure drop vs flow relation can 
lead either to excursive flow or to oscillatory flow. 



Appendix C 


The upper and lower bounds of the integrals 

The integrals given by Eq. IV-94 , IV-101 and IV-111 can be all 

expressed in the form of _ ^ 
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which integrates in 
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where K is a coefficient and m an exponent. For example, the 


integral given by Eq. IV-111 is 
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However, in view of Eq. IV- 34 and IV-28 we have 
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whence we can express Eq. C-3 as 
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By comparing Eq. C-7 with Eq. C-l it can be seen that they are 
of the same form. 

In view of Eq. C-2, the integration of Eq. C-7 yields: 
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which, after some rearrangement, can be expressed in the form of 
Eq. IV-111, thus 



In order to obtain the upper and lower bound of Eq. IV-94, 
IV-101, IV-111 we note that Eq. C-l can be written as 
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where F is the mean value of F given by 
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whence by the mean value thereom 
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which together with Eq. C-10 yields the upper and lower bounds given 
Eq. IV- 95, IV- 105 and IV-116, thus 

x< ' x 


For example, from Eq. C-7 and Eq. C-12 we obtain 



it can be seen that Eq. C-14 can be put in the form of Eq. IV-116. 
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